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I.  INTRODUCTION 


This  first  seal-annual  progress  report  contains  a  senary  of  the  prob¬ 
lem  that  is  being  addressed  in  this  project,  a  senary  of  the  work  that  is 
now  in  progress,  a  list  of  personnel  who  are  participating  in  this  project, 
and  a  senary  of  project-related  activities. 

The  goal  of  this  project  is  to  foranlate  and  investigate  new  approaches 
for  forming  images  of  radar/sonar  targets  from  spotlight-mode,  delay- 
doppler  measurements.  Initially,  we  are  studying  a  particular  processing 
motivated  by  an  approach  used  in  radionuclide  imaging.  Our  longer  term 
goal  is  to  develop  new  processing  based  upon  a  realistic  model  for  the  data 
acquired  with  a  radar-imaging  system. 

Inverse  synthetic-aperture  imaging  (ISAR)  in  radar  and  sonar  relies 
upon  the  relative  motion  between  the  transmitter,  target,  and  receiver.  In 
the  usual  approach,  the  target  is  illuminated  by  a  series  of  transmitted 
pulses.  The  return  for  each  pulse  is  a  superposition  of  reflections  from 
various  locations  on  the  target,  with  each  location  affecting  the  pulse  by 
introducing  both  a  delay  and  doppler  shift.  The  returns  are  processed  to 
produce  an  image  of  the  target.  __ 

The  common  approach  is  to  use  the  same  transmitted-pulse  for  each 
illumination  of  the  target.  Bernfeld  [1]  appears  to  be  the  first  to  intro¬ 
duce  the  idea  for  radar  imaging  of  modifying  the  pulse  shape  on  successive 
illuminations.  Ve  are  using  this  idea  of  Bernfeld's.  Be  also  suggested  an 
approach  for  processing  the  reflected  return  pulses  so  as  to  produce  images 
of  the  target;  his  approach  is  based  on  an  analogy  he  observed  to  the 
equations  governing  the  data  acquired  in  the  x-ray  tomographs  currently 
being  used  for  forming  radiological  images  in  medicine. 

One  deficiency  in  Bernfeld's  novel  spproach  is  that,  for  the  analogy  to 
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x-ray  tomography  to  bi  accurate,  the  ambiguity  function  of  the  transmitted 
pulses  must  be  highly  concentrated  along  lines  in  the  delay-doppler  coordi¬ 
nates  and,  moreover,  must  have  a  constant  amplitude  along  such  lines. 

Thus ,  practical  radar/sonar  pulses  having  ambiguity  functions  with  compli¬ 
cated  sidelobe  structures  and  a  nonuniform  amplitude  along  the  principle 
lobe  are  not  accommodated  very  veil  in  the  concept.  The  purpose  of  the 
initial  phase  of  our  study  is  to  investigate  an  analogy  to  medical  imaging 
vhere  this  restriction  is  relaxed.  This  extension  to  Bernfeld's  idea  may 
permit  improved  images  to  be  formed  for  practical  ambiguity  functions. 

The  analogy  to  medical  imaging  which  we  are  attempting  to  exploit  is 
described  fully  in  our  paper  [2],  which  is  included  herewith  as  Appendix  1. 
To  summarise,  the  transmitted  pulse  has  a  oomplex  envelope  E*  f(t),  where 
E(  is  the  transmitted  energy,  and  the  received  pulse  has  a  complex  envelope 
s(t)  given  by 

s(t)  »  |  b(t-t/2,r)E^2f (v)dt, 

where  b(t,r)  is  a  zero-mean,  complex-valued  Gaussian  process  modeling  a 
diffuse  reflection  interaction  at  the  target;  b(t,r)  is  the  instantaneous 
strength  of  the  reflection  at  time  t  and  two-way  delay  r.  In  our  initial 
study,  we  are  assuming  that  the  scatter  process  is  stationary  temporally 
and  uncorrelated  spatially  (i.e.,  a  VSSUS  model  in  the  terminology  of  Van 
Trees  [3].)  The  power  spectrum  of  b(t,r)  at  a  given  delay  x  is  the  tar¬ 
get's  scattering  function  o(r.f).  Por  two  distinct  delays,  say  and  x^, 
the  processes  Mt.r^)  and  b(t,T2)  are  uncorrelated.  For  our  initial  study, 
the  received  pulse  is  first  processed  by  a  collection  of  bandpass  matched 
filters  and  square-law  envelope  detectors,  as  shown  in  [2,  Figure  1,  see 
Appendix  1];  each  bandpass  matched  filter  is  matched  to  a  doppler-shif ted 


version  of  the  trensaitted  pulse.  This  is  s  different  fora  of  front-end 
processing  froa  the  usnsl  ISAS  processing  where  two-diaensionsl  Fourier 
trsnsforas  ere  used.  Our  aotivations  for  the  use  of  the  BPKF-SLED  ere:  1. 
the  BPKF  hes  known  qualities  for  snppressing  the  effects  of  additive  noise; 
end  2,  the  BPKF-SLED  receiver  erises  in  e  fundeaentel  wey  for  estiaeting 
deley  end  doppler  [3].  Vithout  edditive  noise,  the  expected  velue  of  the 
result  of  this  BPKF-SLED  preprocessing  is  e  function  p(r.f)  of  deley  x  end 
doppler  f  given  by 


p(r.f)  -  Ef  f  <j(t  ' ,  f '  )e(r— t  * .  f-f '  )dt  'df ' , 

—os  -00 

where  e(x.f)  is  the  eabiguity  function  of  the  trensaitted  pulse,  which  is 
the  squared  aagnitude  of  the  coaplex  deley-doppler  correlation  function  [3, 
eqn  (10.18)].  Ve  cell  p(x.f)  the  "deley-doppler  power  function;"  it  is  the 
two-diaensional  convolution  of  the  eabiguity  function  of  the  trensaitted 
pulse  end  the  scattering  function  of  the  target. 

Recognizing  the  effect  on  the  eabiguity  function  of  varying  the  linear 
FH  sweep  (chirp)  rate  of  the  signal  is  iaportent  for  our  work.  The  effect 
is  to  sheer  or  tilt  the  eabiguity  function  in  the  deley-doppler  plane  [3, 
p.  290].  We  denote  this  tilt  by  the  pereaeter  9,  which  depends  on  the 
chirp  rate.  While  it  is  not  necessary  to  do  so,  we  ere  assuaing  in  our 
initial  studies  that  the  pulse  shape  is  changed  along  with  the  the  chirp 
rate  in  such  a  way  that  the  only  variation  of  the  eabiguity  function  is  to 
rotate  it  to  an  angle  6  in  the  deley-doppler  plane;  our  aotivation  for  this 
is  to  aaintain  a  close  analogy  to  radionuclide  inaging.  In  order  to 
include  this  chirp-rate  aodulation  in  our  notation,  we  shall  replace  f(t), 
p(t,f),  s(t).  and  a(x,f)  by  f^(t ) ,  Pg(x,f),  s^ft),  and  ag(r,f)  respective- 
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We  may  now  stete  the  delay-doppler  radar  iaaging  problem  that  we  are 
studying  as  follows.  Estiaate  the  target's  scattering  function  froa  data 
acquired  froa  a  series  of  target  illuminations  for  chirp  rates  resulting  in 
N  tilt  angles  0^,  0^,  ....  8  of  the  aabiguity  function  of  the  transaitted 
pulse.  For  our  initial  study  in  which  the  received  signal  is  preproceased 
by  a  collection  of  BPMF-SLED' s,  the  estiaate  is  to  be  formed  by  using  an 
appropriate  algoritha  on  the  N  delay-doppler  power  functions 

•  CD 

Pq(t,£)  -  EtJ  J  <j(T',f')a0(T-c',f-f')dT'df', 

(for  0  ■  0^,  ....  8  )  appearing  as  the  output  of  the  BPMF-SLED ' s . 

Our  approach  is  based  upon  the  use  of  an  algoritha,  called  the  confi¬ 
dence-weighted  (CW)  algoritha,  that  is  used  in  radionuclide  iaaging  where 
an  expression  siailar  to  that  for  the  delay-doppler  power  functions  p0(?,f) 
is  obtained.  With  reference  to  [4],  in  radionuclide  iaaging,  Pe(r,f) 
corresponds  to  the  intensity  of  detected  coincidences,  at  the  various 
observation  angles,  <r(r,f)  corresponds  to  the  intensity  of  annihilations, 
and  a(r,f)  corresponds  to  the  aeasureaent-error  density,  and  the  goal  is  to 
estiaate  the  intensity  of  annihilations  froa  observations  of  coincidences 
when  the  error  density  is  known  froa  calibration  experiaents.  The  steps 
required  with  this  algoritha  are  detailed  in  [2,  see  Appendix  1]. 

Our  plan  for  studying  the  use  of  the  CW  algoritha  for  radar  iaaging 
consists  of  three  aajor  steps  as  follows: 

1.  establish  a  capability  for  simulating  noise-free  return 

signals  froa  siaple  radar  targets,  siaulating  the  BPMF-SLED 


preprocessing,  simulating  the  CW  algorithm,  simulating  the  j 

i 

conventional  ISAR  algorithm,  and  comparing  results;  j 

I 

2.  establish  a  capability  for  simulating  noisy  return  signals  j 

from  simple  radar  targets,  simulating  the  BPKF-SLED  prepro¬ 
cessing,  simulating  the  CW  algorithm,  simulating  the  conven¬ 
tional  ISAR  algorithm,  and  comparing  results; 

3.  formulate  a  model  for  the  noisy  return  signal  in  an  imaging 
radar  system  and  apply  statistical  estimation  theory  to  de¬ 
rive  a  model-based  algorithm  for  forming  the  radar  image,  and 
compare  this  to  the  results  obtained  with  the  CW  and 
conventional  ISAR  algorithms. 

I 

Our  effort  during  the  first  six  months  has  been  toward  accomplishing  the  j 

first  of  these  tasks. 


II.  SUMMARY  OF  WORK  ACCOMPLISHED 

The  goal  of  our  efforts  during  the  first  six  months  of  this  project  has 
been  to  establish  a  capability  to  simulate  idealized  delay—  doppler  radar 
data  and  to  process  the  simulated  data  using  both  conventional  ISAR  pro¬ 
cessing  and  the  CW  processing  suggested  by  that  used  in  radionuclide  im¬ 
aging,  as  discussed  in  (2,  see  Appendix  1].  The  implementation  is 
"idealized"  in  the  sense  that  no  noise  or  random-ref  lection  effects  are 
included. 

During  this  period,  computer  simulations  of  radar  returns  from  simple 
targets  and  the  processing  of  the  returns  by  radar  image-processing  algo¬ 
rithms  have  been  developed,  and  some  preliminary  tests  have  been  performed. 

Two  targets  have  been  simulated,  a  rotating  rough  disk  and  a  rotating 
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rough  sphere.  The  scattering  functions  for  these  test  objects  are  known 
analytically  and  are  shown  in  [2,  Figures  2  and  5.  see  Appendix  1]. 

Two  types  of  radar  pulse-sequences  have  been  included.  They  are 
illustrated  in  Figure  1.  For  the  first,  a  sequence  of  128  linear  FM 
chirped  pulses,  each  having  a  Gaussian  envelope,  is  transaitted,  with  each 
pulse  having  a  distinct  chirp  rate.  The  duration  and  FM  sweep  rate  of  each 
pulse  were  chosen  so  that  the  aabiguity  functions  associated  with  the 
resulting  pulses  would  all  have  the  sane  shape  in  the  delay-doppler  plane 
but  would  have  their  major  axes  rotated  froa  one  pulse  to  the  next  so  as  to 
cover  all  angles  in  the  range  froa  0  to  180  degrees  uniforaly.  The  paraae- 
ters  for  this  first  pulse  sequence  are  given  in  Table  1. 

For  the  second  series  of  radar  pulses  siaulated,  a  sequence  of  128 
stepped-frequency  bursts  is  transaitted.  with  each  burst  consisting  of  128 
short,  separated  pulses;  each  pulse  in  a  burst  has  a  distinct  frequency, 
but  all  the  bursts  are  identical.  Our  aotivation  for  using  the  step¬ 
ped-frequency  wavefora  is  that  it  is  one  used  at  the  Naval  Ocean  Systeas 
Center  in  San  Diego.  The  paraaeters  for  this  second  pulse  sequence  are 
specified  in  Table  2. 

Two  radar  signal-processing  algorithas  have  been  iapleaented  during 
this  six  aonth  period.  The  first  is  based  on  conventional  ISAR  processing, 
and  the  second  is  based  on  the  processing  used  in  radionuclide  iaaging. 

processing  based  on  conventional  ISAR  techniques-  A  digital  simu¬ 
lation  capability  was  developed  for  generating  ISAR  images.  This 
involved  a  study  of  current  ISAR  processing  techniques  and  the  selec¬ 
tion  of  one  approach  for  iapleaentation.  The  literature  (see,  for 
exaaple,  D.  L.  Mensa  [51)  describes  ISAR  digital  processing  in  terms 
of  two-dim .nsional  discrete  Fourier  transformation  techniques;  this 
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SEQUENCE  1:  GAUSSIAN  ENVELOPE,  LINEAR  FM  CHIRP,  CHIRP-RATE  MODULATION 

^Ireq  (t ) 
env  <  t ) 

->  t  •••  — - — — I  t  •••  — - )  t  •••  - >  t 


A  A- 

* ' V  A.  _  _  _  ^ .  V  i.  -  _  _  i  .  .  .  


PULSE  #: 


CHIRP 

RATE  (MHz/s) 


BAND 

WIDTH  (kHz): 


AMBIGUITY  FNC 
TILT  ANGLE  (0> : 


1 

32 

64 

96 

69.8 

344 

482 

344 

+96 

+80 

0 

-BO 

26 

57.8 

1.8 

57.8 

1.4° 

45° 

90° 

135° 

SEQUENCE  2:  STEPPED-FPEQUENCY  WAVEFORM 


|  HI.. | 

J - )  -j-  a  •  •  J..1.  i - >  t  •  •  • 


t  ...  11JL'-*,  ..>  t  ... 


BURST  #: 


64 


96 


CHIRP 

RATE  <MHz/s>: 


3.4 


3.4 


3.4 


3.  4 


BURST  PARAMETERS: 

PULSE  #:  1 


T 1  =85  ns  — > 

T2  =  577  xs  — > 
T3  =  73  ms  — > 


128 


FIGURE  1.  SEQUENCES  USED  IN  SIMULATIONS 
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TABLE  1:  PARAMETERS  OF  SIMULATIONS  USING  A  GAUSSIAN  PULSE  SEQUENCE 


QUANTITY  SYMBOL  VALUE 


target  distance 

R 

12. 5  nautical  miles 

number  of  pulses 

N 

128 

center  frequency 

fc 

3  GHx 

wavelength 

X 

0.1  m 

speed  of  light 

c 

target  rotation  rate 

m 

down  range  resolution 

td 

cross  range  resolution 
(maj  axis)/ (min  axis) 

rc 

of  ambiguity  function 

e 

target  diameter 

D 

64*rd 

nominal  two-way  delay 

Tr 

2*R/c 

delay  resolution 

rt 

2*rd/c 

doppler  resolution 
angle  major  axis  amb. 
fnc.  makes  with  the 

rf 

2*»*rc/A 

delay  axis,  pulse  i 

0 

( i-1) /an 

constant 

al 

2ne 

constant 

a2 

2rr/e 

parameter 

Pi 

al*sin2(0)  +  a2*cos2 (0) 

parameter 

P2 

(a2-al) *sin(20) 

pulse  std.  deviation 

T 

sqrt  (pl*rt/rf  )/2ji 

pulse  duration  (FWHM) 

FWHM 

2.355*(l+lb IT2) /nT) 
p2/(8nT2) 

FM  sweep  rate 

b 

pulse  bandwidth 

BW 

(l+bT^/T 

pulse  spacing 

T2 

total  illumination  time 

Tt  * 

sum[6*FWHM( i) ]  +  (N-1)T2 

total  aspect  change 

W 

«Tt 

EXAMPLE  USED  IN  SIMULATION: 

tc  *  td  *  0.<  m,  e  »  7, 

w  -  8.88xl0-3 

rad/ s 

target  diameter 

D 

38.4  a 

nominal  two-way  delay 

Tr 

154  ps 

pulse  spacing 

T2 

577  ps 

delay  resolution 

rt 

4  ns 

doppler  resolution 

rf 

0.1066  Hz 

constant 

al 

14n 

constant 

a2 

2n/7 

min  FWHM  pulse  duration 

Tmin 

29.2  ps 

max  FWHM  pulse  duration 

Tmax 

204.5  ps 

min  FM  sweep  rate 

bmin 

0.0  Hz/s 

max  FM  sweep  rate 

bmax 

+/-  287  MHz/s 

max  pulse  bandwidth 

BWmax 

14.7  kHz 

total  illumination  time 

Tt 

316  ms  * 

*  assumes  uo  duty-cycle  limitation 

on  transmitter 

total  aspect  change 

W 

0.161  degrees 

TABLE  2:  PARAMETERS  OF  SIMULATIONS  USING  A  STEPPED  FREQUENCY  SEQUENCE 


QUANTITY  SYMBOL  VALUE 


target  distance  R 

nueber  of  bursts  N 

number  of  pnlses/bnrst  n 
center  frequency  fc 

wavelength  X 

speed  of  light  c 

target  rotation  rate  » 

down  range  resolution  rd 

cross  range  resolution  rc 

target  diameter  D 

bandwidth  B 

burst  duration  T3 

pulse  spacing  T2 

pulse  duration  T1 

nominal  two-way  delay  Tr 

total  illumination  time  Tt 

total  aspect  change  V 


EXAMPLE  USED  IN  SIMULATION: 

rc  »  rd  -  0.6  m.  e  »  7,  «  - 


target  diameter  D 
bandwidth  B 
burst  duration  T3 
pulse  spacing  T2 
pulse  duration  T1 
nominal  two-way  delay  Tr 
total  illumination  time  Tt 
total  aspect  change  f 


12. 5  nautical  miles 
128 
128 
3  GHx 
0.1  m 


64*rd 

o/(2*rd) 

X/  (2*N*rc*w) 

T3/(n-l) 

4*D/c 

2*R/c 

N*T3 

w*Tt 


8.88*10-3  rad/s 

38.4  m 
0.2S  GHz 
73  ms 
577  us 
85  ns 
154  ps 
9.27  s 
4.7  degrees 


i 

.N 

»V 

I 

■  « 
V*. 


i 


n 


V, 

V 


*L 


it  the  approach  that  hat  been  adopted. 

prooeaaint  bated  on  radionnolide  imaalna.  A  digital  simulation 
capability  vat  developed  for  generating  delay-doppler  power  fnnetiona 
for  the  simple  radar  targeta  and  Gauatiait- envelope  radar  pnltea 
detcribed  and  for  prooeating  theae  power  functions  using  the  CV 
algorithm  froa  radionnolide  iaaging. 

Listings  of  the  ooapnter  prograas  that  have  been  developed  to  conduct  these 
siaulations  are  contained  in  Appendix  2. 

Soae  of  our  first  results  obtained  using  the  aiaulation  capability  that 
has  been  developed  are  shown  in  [2.  Figures  2  to  6,  see  Appendix  1].  The 
reconstructions  of  the  scattering  functions  of  the  two  test  targets  ob¬ 
tained  with  our  iapleaentation  of  both  processing  approaches  look  quite 
close  to  the  actual  scattering  functions.  This  has  helped  to  verify  the 
correctness  of  our  iapleaentations.  The  conventional  ISAR  processing  does 
have  sidelobe  artifacts,  with  substantial  power,  that  are  not  present  with 
the  new  approaoh  we  are  studying.  These  nay  be  in  part  due  to  the  fact 
that  we  have  not  yet  introduced  any  windowing  into  the  algorithm  in  order 
to  maintain  the  greatest  resolution.  Windowing  is  now  being  introduced. 

It  is  too  early  in  our  efforts  to  draw  any  major  conclusions,  but  a 
tentative  one  for  the  simulations  we  have  implemented  thus  far  is  that 
conventional  ISAR  processing  and  CW  processing  produce  similar  results  when 
no  additive  noise  or  multiplicative  randomness  sre  present.  We  are  en¬ 
couraged  by  this  beoause  it  indicates  that  CW  processing  can  work.  Our 
expectation  is  that  when  we  introduce  additive  noise  snd  a  random  re¬ 
flection  process,  the  CW  algorithm  will  outperform  the  conventional  ISAR 
processing  both  because  of  the  use  of  the  BPMF-SLED  preprocessing  and  the 
fact  that  the  processing  used  in  rsdionuclide  imaging  was  derived  with 
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Activities  preliminary  to  the  formulation  of  an  optical  ISAR  imaging 
problem  have  begun.  Investigation  to  date  (by  K.  Krause)  is  seeking  to 
determine  the  feasibility  of  extending  the  concepts  employed  at  microwave 
frequencies  to  the  optical  domain.  This  determination  of  feasibility  has 
a  number  of  interesting  conceptual  aspects  that  make  the  problem  distinct. 
For  example,  the  potential  sources  of  image  degradation,  such  as  atmos¬ 
pheric  turbulence,  and  the  particular  effects  seen  in  optical  radar  images, 
such  as  speckle,  need  to  be  considered. 
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III.  VORI  PLANNED 

Additional  work  remains  to  complete  the  first  task  outlined  above  in 
the  Introduction.  Window  functions  need  to  be  introduced  into  the  conven¬ 
tional  ISAR  processing.  The  CW  algorithm  needs  to  be  used  for  the  stepped- 
frequency  waveform,  which  requires  an  evaluation  of  the  requisite  ambiguity 
function.  And  more  comparison  studies  are  needed  to  reach  firm  conclusions 
about  the  relative  performance  of  the  ISAR  and  CW  algorithms.  This  work  is 
in  progress. 

We  have  initiated  an  effort  to  incorporate  additive  noise  and  a  random 
reflection  process  into  our  simulation  of  radar-return  data.  When  this  is 
accomplished,  we  will  process  the  simulated  data  using  both  the  conven¬ 
tional  ISAR  algorithm  and  the  CW  algorithm  in  order  to  compare  results. 

We  have  begun  to  formulate  the  radar-imaging  problem  as  one  of  statis¬ 
tical  estimation.  We  have  begun  to  examine  the  application  of  our  approach 
to  laser-radar  imaging. 
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C.  Lewis,  K.  E.  Krause,  J.  O'Sullivan,  and  J.  T.  Vohlschlaeger.  Soae  brief 
biographical  inforaation  on  each  follows. 


NAME:  Robert  C.  Lewis 

PERSONAL  INFORMATION: 

Birthdate:  June  22,  1955 
Birthplace:  Meaphis,  Tennessee 

EDUCATION : 

BSEE  Meaphis  State  University,  1978 
MS EE  Washington  University,  in  progress 

EMPLOYMENT: 

McDonnell  Douglas  Corporation,  St.  Louis,  1978  to  present 

-  contributed  to  the  research  and  developaent  in  avionics  of 
several  advanced  aerospace  projects  including  the  F-15 
fighter  and  the  Toaahawk  cruise  pissile 

PUBLICATIONS: 

D.  L.  Snyder,  H.  J.  Whitehouse,  J.  T.  Vohlschlaeger,  and  R.  C. 
Lewis,  ' 'A  New  Approach  to  Radar/Sonar  Iaaging, ' '  Proc.  1986  SPIE 
Conf.  on  Advanced  Algorithas  and  Architectures,  Vol.  696,  San 
Diego,  CA,  August  1986. 


NAME:  Kenneth  E.  Krause 

PERSONAL  INFORMATION: 

Birthdate:  Noveaber  3,  1951 
Birthplace:  Decatur,  Illinois 

EDUCATION: 

BSEE  University  of  Illinois,  Urbans,  IL,  1973 
MSEE  Bradley  University,  Peoria,  IL,  1983 
DSc  EE  Washington  University,  in  progress 

additional  coursework  in  Physics,  Math.,  and  Chea. : 

Illinois  State  University,  Noraal,  IL  (30  hours  coapleted) 
1979 


HONORS  AND  AWARDS 
Phi  Kappa  Phi 


EMPLOYMENT: 

General  Telephone  of  Illinois,  June  1973  -  October  1979 

-  special  service  circuit  equipment  and  transmission  level 
specification 

Caterpillar  Tractor  Co.,  Ootober  1979  -  June  1983 

-  electrical  standards,  computer  numerical-control 
troubleshooting 

McDonnell  Aircraft  Co.,  June  1983  -  present 

-  mathematical  modeling  and  computer  analysis  of  radiative 
transfer  processes  in  the  IS  and  visible  bands;  visual 
response  modeling 


:  Joseph  A.  0* Sullivan 

PERSONAL  INFORMATION: 

Birthdate:  January  7,  1960 
Birthplace:  St.  Louis,  Mo. 

EDUCATION: 

BSEE  University  of  Notre  Dame,  South  Bend,  In. ,  1982 

MS  EE  University  of  Notre  Dame,  South  Bend,  In.,  1984 

PhD  EE  University  of  Notre  Dame,  South  Bend,  In. ,  1986 

EMPLOYMENT: 

Aug  1986  -  present:  EE  Dep't,  Washington  University,  Visiting 

Assistant  Professor 

Aug  1985  -  May  1986:  Univ.  of  Notre  Dame,  Senior  Teaching  Fellow 
Jan  1982  -  Aug  1984:  Univ.  of  Notre  Dame,  Research  Assistant 
May  1982  -  Aug  1982:  Abacas  Controls,  Inc.,  New  Jersey, 

Electrical  Engineer 

PROFESSIONAL  AFFILIATIONS  AND  HONORS : 

IEEE,  Eta  Kappa  Nu,  Burns  Fellow 

PUBLICATIONS: 

J.  A.  O'Sullivan  and  M.  K.  Sain,  '*A  Theorem  on  the  Feedback 
Equivalence  of  Nonlinear  Systems  Using  Tensor  Analysis,*'  23rd 
Allerton  Conf.  on  Communication,  Control,  and  Computing,  U.  of 
Illinois,  Urbana,  1985. 

J.  A.  O'Sullivan  and  M.  K.  Sain,  ''Nonlinear  Optimal  Control  with 
Tensors:  Some  Computational  Issues,  "  1985  American  Control  Conf. 
Boston,  MA. 

J.  A.  O'Sullivan  and  M.  K.  Sain,  ''Nonlinear  Feedback  Design: 
Optimal  Responses  by  Tensor  Analysis,"  22nd  Allerton  Conf.  on 
Communication,  Control,  and  Computing,  U.  of  Illinois,  Urbana, 
1984. 


page  13 


.  V  -\ 


NAME: 


a 


i' 

A 


J.  Treat  Wohlschl actor 


PERSONAL  INFORMATION : 

birthdate:  August  12,  1960 

birthplace:  Memphis,  IN 
resident  of:  Mesquite,  TX 


EDUCATION: 

e 

BA 

Math  and  Physics,  Austin  College,  1982 

? 

BSEE 

Washington  University, 

1985 

•  ", 

BSSSE 

Washington  University, 

1985 

MSSSM 

Washington  University, 

1985 

v 

MSEE 

Washington  University, 

1985 

Ph.D.  EE 

Washington  University, 

in  progress 

HONORS  AND  AWARDS: 

v  Class  Valedictorian,  Noth  Mesqnite  H.S. 

s  National  Merit  Scholar 

Member  of  Honorable  Mention  Team,  1985  Math  Competition  in 
Modeling 

jjgl  Sigma  Pi  Sigma,  Tan  Beta  Pi.  Eta  lappa  Nn 

EMPLOYMENT : 

>'  Amoco  Research  Center,  Tnlsa,  OK,  Geophysical  Research  Scientist, 

Sommers  1982-1985 


I 


•*. 


A' 


PUBLICATIONS: 

D.  L.  Snyder,  H.  J.  Vhitehonse,  J.  T.  Wohl schlaeger,  and  R.  C. 
Levis,  ’’A  Nev  Approach  to  Radar/Sonar  Imaging,'*  Proc.  1986  SPIE 
Co iif.  on  Advanced  Algorithms  and  Architectures,  Vol.  696,  San 
Diego,  CA,  August  1986. 


*7 

V 


$ 


* 

V 


I  P 


V.  RELATED  PROJECT  ACTIVITIES 

One  paper  has  been  published  during  the  last  six  months.  A  reprint  is 
included  in  Appendix  1. 

R.  C.  Levis,  M.  I.  Miller,  D.  L.  Snyder,  and  J.  T.  Wohl achlaeger 
visited  the  Naval  Ocean  Systems  Center  in  August  1986  to  discuss  this 
project.  M.  I.  Miller  is  a  member  of  the  Electrical  Engineering  faculty  at 
Washington  University.  The  visit  vas  hosted  by  Mr.  Harper  J.  Whitehouse  of 
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NOSC.  D.  L.  Snyder  gave  a  seminar  detailing  the  ideas  and  goals  of  the 
project  for  interested  NOSC  personnel.  Discussions  vith  several  NOSC 
individuals  interested  in  radar  imaging  took  place  daring  the  day.  A  tour 
of  the  NOSC  radar  imaging  facility  concluded  the  visit. 

(.  E.  Krause,  D.  L.  Snyder,  and  J.  T.  Wohlschlaeger  visited  the  Envi- 
ronmental  Research  Institute  of  Michigan  (ERIM)  in  Ann  Arbor,  MI  on  October 
7,  1986.  The  visit  was  hosted  by  Dr.  James  Fienup  of  ERIM.  D.  L.  Snyder 
gave  a  seminar  detailing  the  ideas  and  goals  of  the  project  for  interested 
ERIM  personnel,  and  discussions  vith  several  individuals  interested  in 
radar  imaging  took  place. 

In  October  1986,  D.  L.  Snyder  gave  a  seminar  about  this  project  in  the 
Department  of  Electrical  Engineering  at  Washington  University  in  St.  Louis. 
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A  NOT  APPBOACS  TO  RADAI/SONAI  DIMING 

Dome  Id  L.  Saydor*.  Harper  J.  Vhitehouee**,  J.  Treat  Vohlaehlaeger*.  tad  Hobart  C.  Lewie* 

•  Vaahiagtoa  Oaivaralty,  St.  Loait,  NO  (3130 
**  Naval  Oaaaa  Syatema  Caatar,  Saa  Ditto.  CA  33133 

Aba tract 

Ve  deacrlbe  aa  aaalogy  betveea  iae|ia|  la  delap— doppler  radar/ eoaar  a ad  poal troa-emiaaioa  tocography. 

Tb  1 a  aaggaeta  aev  prosaaalag  algorlthaa  (or  tbe  radar/ eoaar  imagiag  problao  that  map  permit  improved  vlaaali- 
aatioa  o(  targata  (or  praatiaal  ambiguity  (aaetioaa.  A  racoivar  erehiteetare  eoaaiatiag  o(  a  baadpaaa  matched 
(iltar.  earner a- 1 ev  eavalopa  detector,  aad  apeeialiaad  proeeaaiag  ic  propoaed  to  prodaee  iaagea. 

iHtmwtiti 

lever aa  eyathet io-apertare  imagiag  ia  radar  aad  aoaar  raliea  opoa  the  relative  motioa  betveea  the  traaa- 
altter.  target,  aad  receiver,  la  the  aaaal  approach,  the  target  ia  illamiaated  by  a  aeriea  o(  traaamitted 
palaaa.  The  retara  (or  eaeh  palae  ia  a  aaperpoaitioa  o(  radaetioaa  (rom  varioaa  loeatioaa  oa  the  target, 
vlth  each  looatioa  a((eetiag  the  palae  by  iatrodaeiag  a  ahi(t  ia  delay  aad  doppler.  Ia  thia  way,  the  emaula- 
tive  retara  (or  each  palae  ia  a  complicated  mixtare  o(  rataraa  that  ia  iadaeaead  by  the  ahape  aad  reflective 
propertiaa  of  the  target  aad  ita  motioa  relativa  to  the  traaamitter  aad  raeaivar.  Tbe  raage  aad  doppler 
hiatoriea  of  eaeh  o(  the  redeetad  palace  are  proeaaaed  to  prodaee  a  target* a  image. 

The  common  approach  ia  to  aaa  the  aame  traaamitted-palce  ahape  (ox  aach  lllmaiaatioa  of  the  target. 
Berafeld  (1]  appaara  to  be  the  firat  to  iatrodeee  the  idea  (or  radar  imagiag  of  modifyiag  the  palae  ahape  oa 
aacceaaive  illaaiaatioaa.  He  alee  aaggeated  aa  approach  for  proeeaaiag  the  radeeted  palaaa  ao  aa  to  prodaee 
iaagea  of  the  target;  hia  approach  ia  baaed  oa  aa  aaa logy  he  obaerved  to  the  eqaatioae  goveraiag  the  data  ae- 
qaired  ia  the  x-ray  toaographa  carraatly  baiag  aaed  (or  foraiag  radiologieal  iaagea  ia  aedieiae. 

Oae  dafleiaaey  ia  Berafeld' a  aovel  approach  ia  that,  (or  the  aaalogy  to  x-ray  tomography  to  be  aeeareta. 
tha  aabigeity  (oaetioa  of  the  traaamitted  palaaa  aaat  be  highly  eoaeeatrated  aloag  liaea  la  the  del ay-dopplar 
coordiaataa  aad.  moreover .  aaat  have  a  eoaataat  amplitude  aloag  aaeh  liaea.  Thaa,  pxaetieal  radar/aoaar 
palaaa  haviag  aabigaity  (aaetioaa  vith  eaaplieated  aide lobe  atrnctarea  aad  a  aoaaalfora  aaplitade  aloag  the 
priacipal  lobe  are  aot  aocovmodatad  vary  veil  ia  the  eoaeept.  The  pmrpoee  of  oar  paper  la  to  aote  aa  altaraa- 
tlve  aaalogy  to  medical  imagiag  vhera  thia  reatrietioa  ia  relaxed.  Thia  extaaaioo  to  Berafeld’ a  idea  may 
permit  Improved  imagea  to  be  formed  for  practical  ambigeity  (aaetioaa. 


Tha  model  va  aaa  for  radar/aoaar  imagiag  ia  the  oae  deeerlbed  by  Vaa  Traaa  [3.  Ch.  13].  Tha  traaamitted 
palae  haa  complex  amplitade  K^2f(t),  vhara  ia  tha  traaamitted  eaergy.  Tha  reflected  palae  at  tha  raeaivar 
haa  complex  amplitade  a(t),  chare 


u>  -  jyt 


1  /  % 

1.  *f(t— r)b{t— c/2,r)dt. 


ahera  b(t.t)  la  a  aero-maaa.  complex-valued  Gaoaaiaa  proeeaa  medal  lag  a  diffuae  reflection  iateraetioa  at  tha 
target.  Ia  mill  deaote  tha  power  apeetrm  of  b(t.t)  at  raage  T  by  o(f,t);  thia  ia  the  target* a  acattariag 
faaetioa.  Tea  Trace  deaotee  thia  faaetioa  aa  S^(f,t).  It  ia  eaaamad  that  b(t.Tj)  aad  »(i,tj)  are  uacorra- 

latad  for  A  tj.  Thaa,  b(t,t)  model  a  vida-aeaae  atatioaary,  uacorrelated  aeattariag.  Tha  complex  eavalopa 

of  tha  raealved  algaal  ia  r<t)  ■  aft)  ♦  aft),  vhara  wft)  ia  a  complex-valued,  white  Gauaaiaa  uoiaa-proeeaa 
that  ia  iadapeadeat  of  bft.v)  aad  which  haa  apectral  iateaaity  NQ,  For  the  work  to  be  daaerlbad  hare,  oa 


aaaiaia  that  rft)  ia  firat  proeaaaed  by  a  beady 


tched-filter  eqeare-law  emrelope  deteetor  ( BFNF-SLBD) ,  in 


which  eaae  tha  avaraga  value  of  tha  raaalt  ia  givea  by  [2,  aqa  (13.79)]  pff.t)  ♦  NQ,  where 


pfr.f)  -  B  /P  T*.f*)e(T-t’.f-f’)dV'df 


Thia  work  waa  eupportad  by  the  Office  of  Naval  Beaaarch  aader 
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ia  which  (t.f)  ara  the  delay-doppler  coordiutu  and  a(t.f)  ia  tka  ambiguity  fnaction  of  tka  tranaaittad 
signal,  akick  ia  tka  squared  aa|aitnda  of  tka  coaplez  delty-dopplcr  correlation  faaetioa  [2.  aqn  (10.18)].  Io 
tka  absence  of  additira  noiaa  (i.e.,  ■  0),  tka  average  output  of  tkia  receiver  processing  ia,  therefore, 

tka  tao-diaaaaioaal  convolution  of  tka  target's  scattering  faaotion  with  the  ambiguity  fnaction  of  the 
traneaitted  a 1  goal. 

Adopting  Barnf eld's  idea,  we  eoaaider  that  the  traneaitted  pnlaaa  which  illnainate  the  target  ara  chirp- 
rate  aodalated,  bp  which  wa  aaaa  that  each  pnlaa  haa  a  partienlar  linear  PM  chirp  rata  bat  the  chirp  rate  ia 
wariad  froa  pnlaa  to  pnlaa.  The  effaat  of  changing  the  chirp  rate  of  a  pnlaa  it  to  ahaar  or  tilt  ita 
aabigaity  function  ia  the  delsy-doppler  ooordinatae  (2.  p.  290).  fa  (hall  denote  this  tilt  bp  the  paraaeter 
0.  which  depeada  on  tha  chirp  rata,  fa  will  aasae  ia  what  followt  that  the  pnlaa  thape  ia  changed  along  with 
the  chirp  rate  in  each  a  wap  that  the  only  variation  of  tha  ambiguity  fnaction  ia  to  rotate  it  to  an  angle  0 
ia  tha  dalap-doppler  plana;  it  ia  uaaeee (tary  to  aaiataia  a  fined  thape  of  tha  aabiguity  fnaction  ia  general, 
but  thia  aeaaaptioa  raaaltc  in  tha  oloceat  analogy  to  amission  tonography.  In  order  to  inalade^thia  chirp- 
rate  nodal  at ion  in  on*  notation,  wa  (hall  now  replace  f(t),  a(t.f),  p(t,f),  and  r(t)  by  1^(0.  afl( z.f). 

p#(r.f),  and  r#(t). 

fa  aay  now  atata  tha  delay-doppler  iaaging  problaa  for  radar  and  aoaar  at  that  of  eatiaatiag  tha  target ‘a 

(catteriag  function  frea  tha  data  r_(t)  acquired  for  a  series  of  target  illaainatioaa  at  0  »  6  ,  ....  0  .  By 

"  1  2  N 

teaporarlly  introducing  tha  ciaplifying  aaaaaptioaa  that  the  additive  noiaa  ia  negligible  and  that  the  re- 

aalting  data  can  bo  replaced  by  tkair  average  value,  p^lt.O,  tha  iaagiag  preblaa  baeaaea  a  deconvolution 

problea  ia  which  wo  ara  given  tha  aabiguity  function  a^ft.f)  and  data  pft(t.f)  tor  aovaral  valnaa  of  0,  and  we 

aunt  solve  (1)  for  tha  acattariag  fnaction  a(t,f).  fa  have  noted  ia  CSI  that  thia  iaagiag  problaa  it  anal o- 
goaa  to  one  encountered  ia  poaitron-eaiaaioa  tcaography.  • 

Nodal  for  Potltroa  Halation  TomPtrsphraad  Iaaalna  Problaa 

An  equation  having  the  taae  fora  aa  (1)  oooura  ia  poaitron-aaiatioa  toaographic-iaaging  ayatcaa.  Ia  theta 
ayateaa,  a  positron-emitting  radionuclide  ia  introduced  inaide  a  patient,  and  the  reanltiag  aetivity  ia  ob- 
aerved  eztaraclly  with  aa  array  of  aeiatillatioa  deteetora  that  tnrroand  the  patient  in  a  planar,  ring  geome¬ 
try.  When  a  positron  ia  prodaeed  in  a  radioactive  decay,  it  aanihilatea  with  an  electron  producing  two  311 
kev  photons  that  propagate  at  tke  velocity  of  light  ia  oppocite  direetioaa  along  a  lino.  Ia  ayatana  under 
currant  development,  both  the  liaa-of- flight  and  the  differential  tlme-of-flight  of  the  annihilation  photons 
are  aeaaured,  which  provides  aa  estimate  of  the  location  of  the  annihilation.  Tke  result  la  that,  ia  the 
absence  of  noise,  the  measurements  are  in  tka  fora  of  (1)  with  w(z.f>  being  the  two-dimensional  aetivity 
distribution  to  be  imaged,  and  with  eft(v.f)  being  the  known  error  density  associated  with  the  measuring  the 

location  of  aa  annihilation  event  (4).  for  present  aystaaa,  tha  error  density  is  a  two-dimensional  elliptl- 
s a l- shaped  Gaussian  faaetioa  with  tha  long  axis  oriented  at  an  angle  0  corresponding  to  the  line-of- flight  of 
the  annihilation  photons,  and  data  are  collected  for  from  32  to  96  discrete  eagles  over  the  range  (0,180*). 
Ibis  density  corresponds  to  the  ambiguity  function  of  a  radar  pulse  having  aa  envelope  that  ia  a  Gaussian 
function  sad  a  phase  that  ie  a  linear  ffl  chirp. 

Ia  summary; 

a.  For  delay-doppler  radar/soaar  imaging,  we  suppose  that  the  target  is  illamlnated  by  s  sequence  of  radar 

poises  each  having  a  distinct  PM  chirp  rate  correspoadiag  to  angles  0j,  0^ . 0^  spanning  the  range  from  0 

to  180*  ia  the  delay-doppler  plane.  A  BFNF-SLED  receiver  produces  p0(t,f).  The  smbigulty  fnaction  aQ(T,f)  is 
known.  The  problem  is  to  estimate  the  target's  acattariag  function  o(r.f)  using  the  relationship  (1). 

b.  For  emission-tomography  imaging  when  both  time-of-flight  and  1  lae-of-fl ight  data  are  collected,  we  have  as 
neaewraneata  p((t,f)  for  saglee  0^,  0j,  ....  0^  spaaalag  0  to  180*.  The  measurement-error  density  a#(t,f)  is 

known.  Tha  problem  ie  to  estimate  the  activity  distribution  w(z,f)  using  the  relationship  in  (1). 

Computer  algorithms  for  solving  the  iaaging  problem  for  positron-amission  tomography  with  time-of-flight 
aeasuraaeutt  have  been  under  intensive  development  since  about  1980,  sad  they  continue  to  evolve  and  be  re¬ 
fined.  The  first  such  algorithm,  called  the  "confidence  weighting  algorithm"  {4]  continues  to  be  the  one 
routinely  used.  Algorithms  under  current  development  are  based  on  the  joint  naziaisation  of  likelihood  and 
entropy  (3,61;  these  are  solved  iteratively  and  are,  therefore,  very  demanding  computationally,  bnt  the 
results  appear  to  be  superior.  Ia  the  following  section,  we  indicate  how  the  confidence  weighting  algorithm 
would  be  applied  for  the  radar/soaar  iaagiag  problem.  A  future  publication  will  describe  how  the  newer 
approaches  might  be  used  for  radar  isMglag. 

Coaf Idoaco-Ielahted  Alaorlthn 

The  coaf tdeace—weighted  algorithm  for  procettiag  p#(t.f)  to  estiaate  e(r.t)  is  as  follows.  There  are  two 
steps;  since  the  algorithm  is  linear,  these  could  be  combined  into  one.  but  two  are  used  ia  emission  tcao- 
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inptr  kietui  of  loplMMtttlN  ooiiifwitioai.  The  flrat  (top  1*  to  tors  a  tto-diwuloul  eoavolatioa  of 
tht  data  Pg(t.f)  obtained  for  task  walaa  of  S  with  a  weighting  function  and  thaa  to  foxa  the  pi an 1 -by- 

pi  xnl  ana  of  tkc  reunite  for  task  6;  tkat  la.  wa  fora  tha  foaetioaa 


fe(t.f) 


jj  p#(t\f*>*e(t-r\f-f*)dT'df* 


a  ad  tkaa  wa  oktaia  tka  two-diaanaioaal  ’ 'pr eiaaga ' '  f(r.f)  aecordiag  to 


f(t.f)  •  (1/*) 


J!v- 


f)d». 


Ika  fozaatioa  of  tkia  preinage  general  ires  tka  baek-proj aetioa  atap  of  tka  "(aaf iltarad)  baek-proj aetioa.  poat 
twe-dlaaaaioaal  filtering"  approach  to  Idoalisad  llaa-iatagral  (i.a.,  Badoa  traaafora)  toaography;  for  idaal- 
liao  iatagral  toaography  tha  alaoet  nalweraally  aaad  approach,  which  y ialda  ha  aaaa  flaal  raaalt.  ia  to 
filter  bafora  back  projaotiag,  bat  tkia  doe  a  not  work  wall  whoa  line  iatagrala  arc  a  poor  approxiaatioa.  aa  ia 


tiaa-of-flight  toaography  becaoee  of  tha  diaporaad  poiat  tpraad  pft(x , f) . 


Aa  daaeribod  ia  (31  aad  [4],  wariosa  woigktiag  faactloaa  aight  ba  adopted,  bat  for  tiaw-of-f light  taw 
graphy,  tha  choiaa  aaad  ia  w^(r.f)  ■  a#(r,f).  For  radar/aoaar  iaagiag,  tkia  eorraapoada  to  takiag  the  valao 

of  tka  BFNF-SLED  oatpat  p^U.f)  for  a  partiaalar  walaa  of  (r,f)  aad  diatribatlag  it  owor  tha  delay-doppler 
plaao  according  to  tha  aablgaity  faactioa.  Tkia  aight  ba  aotiwated  by  noting  that  if  wo  ia'tarprat  o(r.f)  aa 


tha  total  reflectance  power  at  (r,f)  aad  a^(r *-r , f *-f)  aa  tka  fraction  of  raflactaaea  power  at  (r,f)  appearing 


aa  aaaaarad  power  at  (t'.f),  tkaa  tha  qaaatity  p  (r.flt'.f)  defined  by 

9 


p#(t.flt*,f*)  -  ae(f-r.f-f)o<t,f)[j^  j^a^r'-t.f'-fJoKr.fldtdf]-1 


ia  tha  fraction  of  raflactaaea  power  at  (r,f)  giwaa  the  total  aaaaarad  power  at  (t*,f),  aad 


o<t.f)  « 


(1/a)  OCL  Pa(r.flT,.f)pe(r,.f')dr'df]d* 


la  tka  total  roflaataaoa  power  at  («,f)  giwaa  tha  power  aaaaarad  at  awary  (t'.f*)  aad  awary  walaa  of  0.  Slaca 
a  property  of  tha  aablgaity  faaetioa  la  that  it  baa  aait  wolaaa  (3,  page  301],  thia  axpraaaioa  ia  aiaply 
another  nay  to  writa  tha  triwial  eqaality  that  o(r,f)  ■  o(r,f).  Bowawer,  if  thia  axpraaaioa  la  adopted  aa  a 
beat  a  for  eetlaatiag  aa  unknown  aeattariag  fa 


action  o(t,f)  giwaa  tha  aaaaaraaaata  pfl(t,f).  and  if  wa  farther 


praaaa  that  tha  aeattariag  faactioa  aaboddad  iaplicitly  within  tha  right-hand  aide  of  thia  axpraaaioa  ia 
eqaalty  likely  to  bo  aay thing  (i.a.,  aalfora)  bafora  aaking  any  aaaaaraaaata,  then  tkia  axpraaaioa  rcdacaa  to 
the  preiaego  axpraaaioa  defined  bafora. 


Tho  aeeoad  proeaaaiag  atap  ia  to  obtain  tha  target'e  iaaga  froa  tha  preiaage.  The  raqairad  operation 
eorraapoada  to  tho  filtering  atap  of  tho  * 'back  project  tkaa  filter'*  aigoritha  aaad  to  iawart  tha  Kadoa 
traaafora.  Tho  analogy  to  tha  coaf ideaee-wel ghting  aigoritha  of  poaitroa-aaiaaloa  toaogrephy  aaggaata  that 
thia  ahomld  ba  aacoapliabad  to  within  a  raaolatioa  faactioa  h(r,f),  which  dafiaaa  a  "daairad  iaaga"  aecordiag 
to 


d(w 


■»-JI 


h(t-x  *,f-f  *)o(r  f '  )dt  'df  ’ 


*a  hawa  toaad  that  iaeladiag  each  a  raaolatioa  inaction  ia  iaportaat  ia  proeaaaiag  aaiaaioa-tcaography  data  aa 
a  way  to  trade  off  raaolatioa  aad  aaootkiag  for  noiaa  aappracaioa.  la  [4],  a  narrow,  two-diaanaioaal ly. 
circularly  ayaaatrlc  Qaaaaiaa  raaolatioa  filter  ia  aaad.  Lot  d(x,f)  denote  tha  aatiaata  of  d(r,f)  obtalaed  by 
proeaaaiag  tho  proiaage  f(w.f).  Alao,  lot  0(a,w)  aad  F(a,w)  denote  tha  two-diaeaaioaal  Fourier  traaaforaa  of 
d(t.f)  aad  f(w.f),  raapaatlwoly.  than,  froa  (4,  eqn .  (17)], 


D(a,w)  -  B(a.w)F(a,w)/0(a,w) , 


nbara  B(a,w)  la  tha  traaafora  of  h(r,f)  aad  0(a,w)  ia  tha  traaafora  of  tha  faactioa  g(r,f)  dofiaad  aecordiag 
to 
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l(t.f)  «  (l/elj^agfT.fJw^tt.flda. 


The  image  d(t.f)  is  obtained  froa  D(u,v)  by  a  two-dimensional  inverse  Fourier  transformation.  Tba  function 
g(t.f)  aad  G(o.t)  ara  prscomputable  aiaea  they  dapoad  only  on  tba  ambiguity  faactloa  aad  tba  weighting 
function  adopted  to  f ora  tba  preimage  aad  not  oa  tba  aaaaarad  data.  For  tba  cbolea  »j(t,f)  •  a^tt.f).  tba 

faactloa  |(t.f)  la  tba  average  orar  0  a f  tba  square  of  tba  ambiguity  faactloa.  Ia  [4],  a^(t,f)  la  a  two- 

dlaaaaioaal  asymmetric  Oaaaalaa  faactloa,  aad  |(t,f)  la  a  Baaaal  faactloa.  Tba  derivation  ia  [4]  doaa  aot 
require  that  a  (t,f)  ba  Gaaaalaa,  bat  thaa  g(t,f)  will  aaaally  need  to  be  evaluated  aaaerically  for  practical 

ambiguity  faaetioaa. 


Tba  arcbitactara  of  the  imagiag  alporitha  wa  bawa  daaeribad  of fare  opportaaitiaa  for  parallel  In  aad  tba 
aaa  of  ipeclally  designed  procacaora  for  aaar  raal-tiaa  processing.  Data  aeqairad  for  aacb  cbirp  rate  9  caa 
ba  procaasad  ia  parallel  aad  then  saasd  to  fora  tba  preimage  f(T.f).  Tba  convolutions  raqairad  for  aacb 
cbirp  rata  aaa  ba  iapiaaaatad  ia  systolic  or  other  high  spaed  architectures.  Tba  architecture  is  dlapraaaad 
ia  Fig.  1. 


fa  bawa  iapiaaaatad  a  preliaiaary  ooaputar  aiaalatioa  to  bagia  tooting  our  algorltha.  Two  targets  bawa 
baaa  slant  a  ted.  a  rotating  rough  dish  aad  a  rotatiag  rough  sphere.  Only  aweraga  quaatitlaa  wars  included  la 
this  first  simulation;  wariations  das  to  both  additive  noise  and  statistical  variations  ia'tha  reflection 
process  ara  prasaatly  baiag  iapiaaaatad.  Processing  was  performed  using  both  oar  new  algorltha  (Cl)  aad  the 
algorltha  conventionally  used  for  inverse  synthetic-aperture  radar  (ISAk)  iaegiag. 

dish:  Tba  aaia  of  rotatioa  is  aoraal  to  tba  disk  aad  passes  through  its  geoaatric  canter.  Tba  radar 
looks  la  tba  plaaa  of  the  disk.  Tba  diaaatar  of  the  disk  is  (.4  asters.  Tba  rotatioa  rata  is  selected  to 
give  a  circular  acattariag  function  ia  tba  delay— dopplar  plane. 

snbara:  Tba  sals  of  rotatioa  it  through  its  gecaatrio  canter.  Tba  radar  looks  ia  the  equatorial  plane 
of  the  sphere.  The  diameter  of  tba  sphere  it  4.4  aetara.  The  rotatioa  rats  is  selected  to  give  a  circular 
acattariag  function  ia  tba  delay-doppler  plaaa. 

Figures  1(a)  aad  Kb)  show  tbs  power  acattariag  function  for  the  rotatiag  disk  and  its  reconstruction 
using  the  Cf  algorltha.  respectively .  Each  laags  ia  aa  array  of  128  range  by  US  urosa-raaga  resolution 
calls,  with  each  cell  baing  a  square  of  dlacnaion  0.1  aatar  oa  a  side.  Tba  scattering  function  is  coatelaad 
in  a  square  aabarray  of  <4-by-<4  piasls.  A  total  of  12S  distinct  FM  chirp  rates  vase  used  such  that  the 
eabiguity  function  rotated  at  equal  incraBsats  through  180*.  Tba  asymmetric  Gaussian- shaped  ambiguity 
faaetion  bee  dimensions  of  28  piaals  (FWHM)  along  tbs  aajor  axis  aad  4  pixels  (FVHM)  along  tba  minor  axis. 

The  durations  sad  chirp  rates  of  the  radar  pulsss  wars  adjusted  so  that  tbs  allipsas  describing  their  ambigui¬ 
ty  fuactioaa  had  tba  same  aajor  end  minor  axis-diaanaioss  at  each  rotation  sugla.  The  delay-doppler  iasge 
obtained  with  tba  Cl  algorithm  shows  same  degradation  ia  resolution  but  no  other  major  distortions  are  evi¬ 
dent. 

Flguree  3(e)  aad  3(b)  shoo  the  power  scattering  function  for  the  rotating  sphere  and  its  reconstruction 
nsiag  the  Of  algorithm,  respectively.  Each  image  is  sa  array  of  128  rsnge  by  128  cross-range  pixels,  with 
each  pixel  baiag  a  square  of  dimeuaiou  0.1  mater  oa  a  side.  Tba  scattering  function  is  contained  in  a  square 
subarray  of  44-by-44  pixels.  A  total  of  128  distinct  FM  chirp  rates  were  used  such  that  tba  ambiguity 
faaetioa  rotated  through  180*.  The  durations  aad  chirp  rates  of  tba  radar  paisas  ware  adjuated  so  that  the 

allipsea  deearlbiag  thair  ambiguity  functions  bed  the  same  major  and  minor  axis-dimensions  (24  by  4  pixels 

FVIM,  respectively)  at  each  rotatioa  angle.  Tba  delay-doppler  lauge  obtained  with  the  Of  algorltha  shows  some 
degradation  ia  resolution  but  ao  other  aajor  distortions  are  evident. 

Figurae  4(e),  (b),  aad  (c)  show  examples  of  pQ(r,f)  for  9  •  0* ,  43*.  end  90*  for  the  rotating  sphere. 

Figuree  3(e)  end  3(b)  show  the  squared  megaitude  of  the  reflectivity  (complex  amplitude)  function  for  the 
rotatiag  disk  and  ite  reconstruction  using  the  conventional  ISA!  algorithm.  Each  image  is  sn  srrsy  of  128 
reaga  by  128  eroes-raage  resolatioa  calls,  with  each  call  baiag  a  square  of  dimension  0.1  metar  oa  a  side. 

Tba  reflectivity  faaetioa  ie  coataiasd  ia  a  square  subsrray  of  32-by-32  pixels.  A  stappad-freqsancy  radar- 
pulse  eequaaca  eoaeistiug  of  128  square  pulses,  each  of  duration  83  as  sspareted  ia  time  and  frequency  by  3.4 
me  sad  11. T  Hi.  raspastivsly.  is  assmsd.  the  center  frequency  is  3  GBi,  sad  tbs  target  is  illuminated  by 

128  such  pulse  sequences.  For  tbs  simulation,  we  calculated  tbs  received  signal  hem  the  target  using  dis¬ 

crete  points  to  model  the  distributed  reflectivity  within  resolution  cells  distributed  evenly  ia  a  reetengnlar 
grid  over  the  target.  The  received  signal  for  each  traaMitted  pulse  is  processed  to  generate  the  image  by 
first  sorting  into  range-resolution  cells  using  a  Fourier  transform  and  then  sorting  the  result  into  doppler 
resolution-cells  by  again  aaiag  a  Fourier  transform  oa  tbs  data  in  each  range  cell.  The  resulting  image  shows 
some  loss  is  resolatioa  sad  aa  oscillatory  artifact  around  the  edge  of  the  disk. 
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Flgarta  6(»)  aad  6(b)  «ho»  tka  u|ail(4i  of  tkt  ti<Ucli«itr  (uitloa  (or  tbo  tatitlii  i|km  itk  Its 
racoat true t ion  axing  tbo  cauattieul  lilt  algoritkaf  bck  int|i  in  nn  array  of  12g~by-124  rnnointion  cnlln. 
Tit  rtflactivity  function  it  eontainod  in  t  (quart  anbarray  at  32-by-32  pixala.  »itk  tack  pixal  baing  a  tquart 
of  diaaaaioa  0.1  anttr  on  a  aid*,  lit  ta at  ISAM  proaoaaia |  daaeribod  for  tka  diak  van  trad.  Tka  ratal  tiny 
inaga  akoaa  lota  in  ratolatioa  and  as  oacillatory  artifaat  aronad  tka  adga  of  tka  diak. 

■a  viafc  to  aapkaaiaa  tkat  tkaaa  raaalta  art  oar  firat  oaaa,  and  tkay  art  vary  pralialaary.  Va  art 
aacoaragad  by  tkaa  baeaaaa  tka  mage  a  prodacad  t itH  tka  Of  algorithm  da  ataa  to  kart  lata  artifact!  tkaa  tkoaa 
prodacad  aitk  tka  aoavaatioaal  ISAM  algoritka.  fa  tapoat  tkat  a  grot  tar  diffaraaaa  la  tka  retaltiag  ioagae 
•ill  raaalt  wkaa  additive  aoi at  and  atatlatieal  variatlaaa  ia  tka  raflaatioa  proetaa  art  iaeladad  in  oar 
aiaalatioa.  aitk  tka  Of  algorithm  baiag  t upt r  1  or  baaaaat  of  tka  aaa  of  tka  aatakad  flltar  for  additive  aoiaa 
lappraaaioa  and  averaging  o»tr  dopplar  rataa  ia  foraiag  tka  iaaga. 
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Figaro  2.  a)  pooor  aeattoring  faaotloa  for  a  ro¬ 
ta  ting  diak.  b)  racoattrnctloa  aaiag  tko  Of 
algorithm. 


Figaro  1.  a)  goaarata  tka  dolay-dopplar  faaetioa 
far  tack  ekirp  rata,  b)  gaaorate  praiaagt  by 
•tigkttd  btck-proj tetloa  aad  tka  dotirtd  aatlaata 
by  poat  filttriag. 
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Figara  3.  a)  power  aeatteriag  faaetioa  for  <  t-~ 
tatfag  iphere.  b)  raeoaetrootioo  aaiag  tb»  a 
•lgoritba. 
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Figaro  4.  BPMF-SLED  dalay-dopplar  oatpot  for 
•  )  0*.  b)  43*.  aid  e)  90*. 
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Figaro  S.  a)  aq  oared  aagaitada  of  tha  reflectivity 
faaetioa  for  a  rotatiag  dial.  b)  rocoaatraetioa 
aaiag  tha  ataadard  ISAS  algorithm. 


Figaro  6.  a)  aqaarod  aagaitada  of  tha  reflectivity 
faaetioa  for  a  rotatiag  aphero.  b)  rocoaatroetion 
aaiag  tha  ataadard  ISAS  algoritba. 
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APPEND LI  2A.  PROGRAM  LISTINGS  POE  CONVENTIONAL  ISAR  SIMULATION 


SIMULATION  PERFORMED  BY:  Robert  C.  Levis 
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A  FORTRAN  prograa  was  written  to  achieve  the  goal  of  prodnoing  conven¬ 
tional  ISAR  iaagea  in  a  laboratory  siaulation.  The  diaenaiona  of  aianlated 
targeta  were  apeoified  and  their  aoattering  fnactiona  were  anpplied  for 
nee.  The  ISAR  aianlation  contained  two  parta.  The  firat  inpnta  the  scat- 
tering  function  and  calcnlatea  the  received  signals  froa  that  target  for 
all  pulses  and  bursts  tranaaitted.  The  second  part  processes  the  received 
signals  using  Fourier  transfora  techniques,  as  described  here,  to  produce 
an  iaage  data  file.  An  iaage  display  routine  is  used  to  display  the  ISAR 
iaage  on  a  color  aonitor. 

An  algoritha  is  presented  in  Figure  A2.1  to  describe  the  prograaaing  of 
the  aizer  output  portion  of  the  aianlation.  A  two-diaensional  array  is 
used  to  store  the  resulting  radar  signals.  Two  analytic  targets  were 
chosen,  a  rough  disk  and  a  rough  sphere.  The  field  of  view  of  the  iaage 
was  chosen  to  be  76. S  aeters  in  both  range  and  oross  range  diaensions.  The 
resolution  in  both  diaensions  was  chosen  to  be  60  ca.  Therefore,  the 
scattering  function  is  a  128-by-128  array  of  reflected  power  values.  The 
diaaeters  of  both  the  disk  and  sphere  were  chosen  to  be  38.4  a  with  both 
centered  in  the  field  of  view.  It  was  desired  to  write  a  coaputer  prograa 
that  would  generate  an  iaage  of  size  256-by-256  pixels.  This  is  required 
by  the  Fourier  transfora  techniques  used. 

In  this  simulation,  128  pulses  per  burst  and  128  bursts  are  used  in  the 
transaitted  wavefora.  This  covers  the  field  of  view  given  the  resolution. 
The  resulting  two-diaensional  array  holding  the  received  radar  signal  is  of 
size  128-by-128.  The  range  profiles  are  calculated  by  using  a  discrete 
Fourier  transfora  (DFT)  and  the  following  property  to  get  an  inverse  DFT: 
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Mila 


R**d  in  Scattering  Function 
Increment  Burst  Number  a 
K*1 

Increment  Resolution  Cell  (r,  f)  •  (i,  j) 
Increment  Pulse  Number  n 
I  Calculate  Received  Signal  S(n) 

I  S.(n)  ■  S.  ,(n)  ♦  S(n) 

End  LBop  ** 1 
K  •  K  ♦  1 
End  Loop 

Store  S(n)  in  S(n.  a) 

End  Loop 

Normalize  S(n,  a)  to  mx  [S(n,  «)]  •  1.0 
End 


FIGURE  A2.1  Algorithm  Used  for  Received  Signal  Simulation 

if  X(k)  is  the  DFT  of  x(n).  then  Nx*<k)  is  the  DPT  of  X*(n),  where 
denotes  complex  conjugation.  Alep,  in  computing  a  range  profile,  the  128 
point  signal  is  placed  in  the  center  of  a  236  point  array,  padded  with 
zeros,  for  the  inverse  DFT  operation.  The  resulting  256  point  range  pro¬ 
file,  divided  by  the  number  of  points,  replaces  the  row  of  signal  data  in 
the  two-dimensional  array. 

In  the  final  step  of  calculations,  the  cross-range  profiles  of  the 
image  are  made  using  a  DFT  on  the  columns  of  the  two-dimensional  array 
resulting  from  the  last  operation.  However,  a  rearrangement  of  the  data  in 
the  236  point  array  is  needed.  Instead  of  placing  the  128  data  points  in 
the  center  of  the  256  point  array,  the  first  64  data  points  are  placed  in 
the  first  64  array  locations,  leaving  the  center  array  locations  zeroed. 


This  is  illustrated  in  Figure  A2.2.  After  the  DFT  is  computed,  the  zero 
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frequency  location  is  in  the  first  array  location,  so  the  data  are  rear¬ 
ranged  with  the  last  half  of  the  256  points  interchanged  with  the  first 
half.  The  magnitudes  of  the  resulting  complex  numbers  are  stored  in  each 
colnsu  of  the  two-dimensional  array.  The  resulting  two-dimensional  image 
data  array  is  then  used  with  a  color  scale  to  display  the  target's  image. 


FIGURE  A2.2  Fourier  Transform  Techniques 


A  listing  of  the  computer  program  follows.  In  the  main  routine,  PR0C3, 
the  target's  scattering  function  is  input  and  the  radar  received  signal  is 
calculated.  This  signal  is  normalized  to  peak  magnitude  in  loops  110,  120, 
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130.  tad  140.  Subroutine  IMASB2D  is  thta  called  to  reconstruct  the  image 
A  brief  description  of  eaob  program  and  sabrontine  is  as  follows. 
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PK0C3: 


oalenlates  a  sixer  output  signal  fraa  an  inpnt  scattering 
fnaetion  data  file.  A  aoarstatistioal  target  is  assmed. 


CSX.  (XT 


TT 

BW 

F0 

0 

N 


location  of  center  of  rotation  within  the 

scattering  fnaetion  in  oa. 

total  tiae  interval  of  all  palsee  and  bursts 

bandwidth  of  a  burst 

lowest  frequency  in  a  burst 

step  frequency  inereaeat 

amber  of  pulses  in  a  burst,  amber  of  bursts 


IMA6K2D:  performs  I&AK  processing  to  reconstruct  an  image  from  radar 
signals.  The  DFTs  are  implemented  using  the  FFT  algorithm. 

TABfiBT:  reads  in  the  target's  scattering  fmetioa  data. 

PROGRAM  PR.0C7 

D I MENS  I  ON  3< 123 > ,9(123) ,31(1 2S ? 

PEAL  LAMBDA , NA 
Complex  FRAME 

Common  /b 1 k 1 /  LAMBDA (123,123) 

C omm on  / b  1  k  2,- "  F RAM E  <.  25 6 , 25 6  > 


C*  *****************************************************  ******* 
C*  STEPPED-FREQUENCY  I SAR  PROCESSING  ( NINCCNEC)  * 

C*  NON-STAT I  ST  I  CAL  REFLECT  I’.'  1 7Y 

C*  BOB  LEU.' I  S  ,  NOV.  1986 

'  *  ;-.put=16,  output=17 

C=3 . E3 

FI =3. 1415927 
C  RX=  6 . 3 

_  z  o 

(7.  I  —  'Z  *  W 

77=9 . 23 
BU=0 .  2CE9 
F9=2 . 375E9 
DF=BW/ 127. 

‘■1=1  23 

■2+ii.*r*:Wit:*>tr#W  +  *-+*:**it:#*******r**-4:**:*#:*******-************:*****~ 

<2********  Read  in  Target's  scattering  -function.  ♦*•*****♦*♦ 

g  ^-s******  ■■t-****#****-*****************-*- 

CALL  target 
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DG  50  K=1 ,123 


S  <  K  >  =0 . 

Q( K ) =  0  . 

SI ( K>  =  <  K-l >  *T2+T +T 1 2 
50  CONTINUE 

DO  80  1=1 ,123 

X= ( I  - 1 >  * . 1 -CRX 
DO  70  J=1  ,  1  23 

I r  ( LAMBDA  < I , J  >  . EQ .3.) GO  TO  73 
A= LAMBDA ( I , J>**9 .5 
Y= < J- 1 ) * . 1 
DELAY=6.6667E-9*Y 
DO  60  K=1 ,123 

FI  =  <  K-l )*DF+F3 
NT2C=X/ (12. 3*DT*FI  ) 

PI 2FTAU=6 . 2332*FI  *<  DELAY -VT2C* 

S(K>=S<K>+A*C0S<PI2FTAU> 

GK  K)=Q(  fO  +A*SIN< PI 2FTAU) 

£3  CONTINUE 

70  CONTINUE 

30  CONTINUE 

DO  90  K- 1,1 28 

FRAME <  L , K+  64 > =CMPLX  < 3  < K  > , Q  < K  > ) 

^0  CONTINUE 

L=L+  1 
130  CONTINUE 


?K=8  . 


DC  120  L* 1,256 
DO  110  K= 1,256 

P K=MAX( PK ,ABS<  FRAME < L , K) > > 
110  CONTINUE 

120  CONTINUE 

K I  NM=  1  ./PK 
DO  140  L* 1,256 
DO  130  K= 1,256 

FRAME  <  L , K  > = FRAME ' L , K > *PK I NU 
133  CONTINUE 


3 b r  *  t  n  ■?  image  2  D 


Real  image 

Complex  + < 256 > , i r ame 
Dimension  UJ1C123) 

C ommon  /b 1 k 2/  f r ame < 256 , 256  > 
Common  /b 1 k 3/  i mage < 256 , 256 > 

M=3 

3F=5.961 
DO  10  1=1,123 

W 1  (  I  )  =3F*EXP  <  -3  .  0 0 I  -64  >  **2  ' 
CONTINUE 


DO  120  i=l ,256 
DO  1Q0  j=l , 256 

•f  <  j  )  =c  on  j  g  ( t  r  ame  <  i  ,  j  '>  ) 

00  CONTINUE 

DO  101  j=l ,123 

*  <  j  ♦  64 .)  =-f  *:  j  +  64  >  «W  1  <r  j  ) 

101  CONTINUE 

CALL  FFT ■'  i  ,M> 

DO  1 10  j=l ,256 

■frame  <  i  ,  j  )=con  j  g(  -f  ■'  j  )  >  *0 .00  390 
110  CONTINUE 

[20  CONTINUE 


i  25 


1  30 


i  40 
150 


DO  150  j  =  1  ,256 
DO  125  i= 1,256 

f  (.  i  )  =cmp  1  x  (  0  .  ,  0  .  > 

CONTINUE 
DO  130  i = 1 ,64 

■f  (  i  +  64 ) -4  r  ame  (  i  + 1  23  ,  j  > 
f  (  i  +  1  23  > =ir-  ame  {.  \  *  64  ,  j  > 
CONTINUE 
DO  131  i = 1 , 1 28 

■f  <  i  +  64 ) =W 1  C  i  >  *  +  <,  i  +  64 ) 

CONTINUE 
CALL  FFT < f ,M> 

DO  140  .=1,123 

i  mage  (.  \  *  1  23  ,  j  'i=  abs<  f  ‘  i  >  >  >  **2 
i  mage  (.  i  ,  .j  )  =  (  abs  (i<  •  *122''  >  '•  *+2 
CONTINUE 
C CNT I NUE 


,  1 


MAGE  < 


\  —  - 


c  a  *  u  r  n 
1  ~ C PM AT ' 1 


3E9. ;■ 


page  29 


SUBROUTINE  TARGET 


C********  READ  IN  TARGET'S  SCATTERING  FUNCTION.  ******** 


REAL  LAMBDA 

COMMON  / 3LK1 /  LAM8D A ( 1 2 8/ 1 28 > 


R£ADC16,1) C (LAMBDA 1*1,1 28), 28) 


RETURN 

1  FORMAT (1X,8E9.2) 
END 


appUMPIT  IB.  FK06KAM  LISTINGS  DOS  COHFMBHCB  WEIGHTED  ISA1  SMJLAXION 


gnim -ATT ON  POTBOMIHP  Bfi  J. 


Treat  Tokl*okla«c*r 


The  software  system  developed  to  simulate  the  proposed 
processing  algorithm  consists  of  several  separate  programs, 
each  designed  to  perform  a  specific  step  of  the  processing. 

The  programs  are  listed  below  in  the  order  that  they  are 
run,  along  with  a  short  comment  on  their  function.  For  a  more 
thorough  desription  of  a  particular  program's  use,  please  see 
the  accompanying  listings. 


Pr  egr  am 
1 .  mkf i 1 t 

2  .  s  p  i  n  b  a  1  ’ 


Descr i p  t i on 

Generates  ambiguity  functions.  Run  once  f : "■  a 
given  set  of  parameters. 

Generates  the  scattering  function  for  a  rotating 
sphere. 


3a.  detrtrn  Generates  the  p(theta)'s  corresponding  tc  radar 

returns  from  a  deterministic  reflection  model. 


?  t  .  r  r  dm  r  t  r  n 


4  . 


C 1  V  a  S  d  2  d 


Generates  the  p'  theta)  's  correspond:  r.  g  tc  r  a: 
returns  irom  a  random  reflection  model  . 

Performs  the  confidence  weighted  algorithm:  c 
the  returns  data  generated  in  step  3. 


•Generates  the  final  image  by  filtering  f*e 
p  r  mage  de>  'el  oped  '  n  step  4. 

d  ;  s  p  ’  a  ■/  s  the  final  image. 
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u;»  i  l 


mkf  i  1  t  .  f 


o  gr  am 


:  a *e  Ur  i  tten 
UJr  i  1 1 er  f  or 


J.  "rent  Woh 1 schal eger 

March  21,  1985 

Delay  Doppler  Radar  Imaging 


am  Intent 

This  program  has  several  purposes.  The  miain 
■i  is  to  generate  a  lookup  table  that  will  be  used 
an o t h e r  p r ogr am  to  per  f orm  convolutions.  These 
•  v o I  u  t  i  ons  w  i  1  1  be  w  i  th  a  spec  i  -f  i  c  type  o-f  f  unc  t  i  on  , 
ie 1 >  a  two-dimensional  assyme tr i c  gaussian.  The 
•cial  nature  o-f  this  -function  makes  the  task  o-f 
iding  the  look  up  table  coefficients  a  little  easier 
■faster)  than  -for  an  arbitrary  -function. 

'f  performing  the  convolution,  only  pixels 
>. t  -fall  within  a  certain  range  o-f  the  image  pixel 
-rently  being  convolved  will  be  updated.  This 
ge  is  chosen  to  be  a  certain  number  o-f  standard 

■  at  ons,  where  a.  standard  deviation  is  related 
the  -full  width  half  max  me  asu  remen  t  error  by 

:  or  start.  The  -first  step  is  therefore  to  find 
o-  pixels  fall  within  a  given  number  of  standard 
; at  ions  of  a  central  point.  This  is  done  by 
'.■erasing  an  ellipse  of  isoprobability  with  the 
:ropriate  semimajor  and  semiminor  axes  and  testing 
see  whether  the  center  of  a  pixel  lies  inside 
-ellipse.  If  it  does,  the  pixel's  row  and  column 

■  tabulated.  This  is  accomplished  for  all  "nang" 


numei 
a  b  o ' 1  • 


Then  for  all  angles  fromi  the  first  (at  angle 
'  >  to  the  angle  corresponding  to  theta0  +  45 
■s,  for  all  of  the  pixel  centers  inside  of 
1 i p se  and  tabulated  earlier,  a  two-dimensional 
cal  integration  is  performed  to  find  the  volume 
the  pixel  under  the  gaussian  function.  The 
al  method  is  a  simple  riemann  surri  with  "numint" 
on  evaluations  in  each  direction  for  a  total 
■  ■  r,  t**2  function  evaluations.  The  volume 
esth  the  gaussian  for  a  pixel  becomes  the 
:  ■ ng  coefficient  or  look  up  table  coefficient  for 
e’ .  Now,  each  pixel  which  is  assigned  a 
cent  must  also  have  an  address  so  that  the 
■rr  which  is  performing  the  subsequent  convolutions 
now  >where  to  apply  this  weight.  The  address  is 
s  ted  by  taking  n  (the  image  array  size)  time  s 
w  r  which  the  pixel  is  located  plus  the  column. 


offset  =  ( n  *  r ow )  +  c  c 1  . 


a  positive  c-  negat 
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*  * 

c 

a  (-2 

,"2) 

a( -1 , -2)  a ( 0 , - 2  > 

a< 1 , -2) 

a(2,- 

2) 

c 

a  <  -  2 

»-:> 

a  ( - 1  ,  - 1 )  a  ( @  - 1  > 

a<  1,-1) 

a  <  2 ,  - 

1 ) 

1 

,  0) 

a < - 1 ,  0 )  a( 0 ,  0) 

a(l,  0) 

a(2, 

0) 

- 

-2 

4 

?  *  ‘ 

a ( - 1 ,  1)  a  <  0  ,  1) 

a( 1  ,  1 ) 

a(  2 , 

1  ) 

V 

c 

c 

a  (-2 

j  2) 

a  (  - 1  ,  2 )  a  <  0  ,  2 ) 

a< 1 ,  2) 

a<  2 , 

2) 

3 

r 

Note 

t  h  a  t 

the  column  index  is 

shown 

before 

the  row 

y 

z 

r-  -J  _ 

•>  ,  *nd 

that  the  column  number  increases 

from 

z 

J.  U( 

left  t 

o  the  right  in  the 

image  and  the 

row  number 

*N 

- 

i  n  c  r 

e  ase  s 

from  the  top  to  the 

bottom 

.  The 

gauss i an 

function  is  assumed  to  be  centered  at  a  (0,0).  The 
address  (or  offset)  of  pixel  a<-l,-2)  would  be 
(n*row)  +  col  =  <5*<-2))  +  <-l)  =  -11.  The 
addresses  of  the  other  pixels  are  defined  similarly. 

Note  that  earlier  it  was  stated  that  the 
numerical  integration  was  performed  only  for  the 
gauss i an  function  generated  for  angles  from  the 
first  (at  angle  theta0>  to  the  angle  corresponding 
to  theta©  plus  45  degrees.  This  is  not  quite 
true.  While  the  integration  is  performed  for  the 
gauss i an  generated  at  only  these  angles,  the  integration 
will  produce  results  that  have  some  symmetry.  The 
volume  over  pixel  a(l,l>,  for  example,  will  be  the 
same  as  the  volume  over  pixel  a<-l,-l>  for  any  given 
angle.  This  program  takes  advantage  of  this  symmetry 
and  performs  the  integration  only  for  roughly  half  of 
the  pixels  falling  inside  of  the  ellipse.  Note  that 
the  number  of  pixels  inside  of  the  ellipse  will  always 
be  odd.  Therefore  < < num i ns< angl e )  +  1 >/2>  volumes 
must  be  taken,  where  num ins  is  the  number  of  pixels 
falling  inside  of  the  ellipse. 

Fine.  Now  we  have  the  addresses  and  coefficients 
for  the  angles  from  the  first  to  the  one  corresponding 
to  the  first  plus  45  degrees.  Now  we  take  advantage 
of  the  fact  that  the  next  45  degrees  worth  of  coefficients 
and  addresses  can  be  found  merely  by  exchanging  the  roles 
of  the  rows  and  columns  found  for  the  first  batch  of 
erg’e  =  .  Similarly,  the  last  96  degrees  worth  of  addresses 
and  coeff  cients  can  be  derived  from  the  first  90  degrees' 
results  by  merely  changing  the  sign  of  the  row.  All  of 
these  symmetry  properties  are  exploited  to  minimize  the 
rur.  time.  Note  that  the  bulk  of  the  run  time  is 
spent  doing  the  numerical  integration,  but  this  need 
be  done  only  for  approximately  1/8  th  of  the  pixels 
for  which  addresses  and  coefficients  are  calculated. 

At  the  beginning  of  the  program  description, 

I  mentioned  that  the  program  has  several  purposes, 
but  I  have  mentioned  only  the  generation  of  the 
c ok  u p  table  ■  t  se 1 f  .  Other  tunc  t i on s  p e r f orme d 
-■’_de  r~.zi t  *  ~  the  =ze  of  the  b=k  file  the 


tab; e  w 


n  o  i  c  a.  t  i  n  g  if  any  c  * 


ts  are  neoati 


( they  shou Id  ell 
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;  d  e  a 1  1  :•  ,  and  ou  tpu  t  *  i  ng  the  number  o-f  pixels  included 

r. s i de  the  ellipse  at  each  angle  and  the  total  volume 
over  these  p  e  '  s  . 


nput  Fa-  arr.e  t  er  s  : 

fwhme  =  full  width  half  max  measurement  error 
•  along  the  line  o-f  -flight,  in  cm. 

twhiTit  =  -full  width  hal-f  max  measurement  error 

transverse  to  the  line  o-f  -flight,  in  cm . 
p i x  =  i mage  pixel  size,  in  cm . 
nostde  =  number  o-f  standard  deviations  worth  o-f 

pixels  -for  which  addresses  and  coefficients 
are  to  be  generated. 

nang  =  number  of  angles  for  which  the  table  is 
to  be  generated. 

nurr, ;  n  t  =  number  of  intervals  in  which  the  gaussian 
function  is  to  be  subdivided  for  the 
numerical  integration,  in  each  direction, 
n  =  image  array  size,  i.e.  the  image  is  of 

size  n  x  n . 

‘ha  first  four  parameters  are  real  floating  point 
: umbers,  and  the  last  three  are  integer  numbers. 


n  a  n  g  = 
n  urn  ;  n  t 


n  p  u  t  files 

u  4  =  parameter  table  (see  above) 


output  files 

lu3  =  look  up  table  containing  addresses  and  coefficients 
of  the  two-dimensi ona?  gaussian  for  all  "nang" 
angles 


a 


program,  structure 
mkf  lit,  newf i o 


3 


3 


res’  c oef f 1 <  1 060 )  ,  coef f 2( 1 000 , 1 28) ,  confac,  cosqth,  cos the , 

*  del,  fwhmb,  f whme ,  nostde,  pi,  pix,  scale,  semaj , 

*  sen-ii. semin,  semisq.  sigbsq,  sigesq,  sigmab,  s  i  gmae , 

*  sin  cos,  sinthe,  sisqth,  theta,  u,  volpix,  volume (128),  x, 

*  xcert,  xlolim,  xmax ,  xm i n ,  xrot ,  xrotsq ,  xsqfac,  xsqmul, 

*  x'/f  a:,  xymu1,  y ,  ycent,  ylolim,  ymaxsq,  ymiinsq,  yrot, 

*  :.<t  o  t  s  q  ,  >  =  q  f  s.c  ,  ysqmu  1 

■*teger*4  angnum,  angl ,  ang2,  col,  colmax,  colmin, 


*■ 

-ol tab 

(  10  0  0  , 

128),  in 

dex ,  i owr  t ,  1 u3, 

n , 

A 

■ 

♦ 

r,  urn, .  n  s 

<128)  , 

num i n  t  , 

numneg , 

numrec  , 

offset(1008) 

r an  add 

,  r  oi'j , 

r  owmax  , 

r  owm i n  , 

r  ow t  ab ( 1 0 6 0 , 1 28 ) , 

*■ 

t  c  t  n  urn. 

,  x  i  n  d 

x  ,  y  i  n  dx 

>• 

• 

r  t e ge r ♦4 

c  u  t  d  e 

s,  nbyte 

s  ,  n  b  1 

j 

n  te  ger*4 

n  an  g , 

n  terms , 

tot  an g , 

m  i  nor 

i 

! 

r,  t  e  g  e  r 

1 1  o  u  t 

,  1 1 :  n 

e  i  ' 

fwhme 

,  -Kh-.hmt  , 

cmatrx ( 

__  1  ho  t  4  0*7 
I  tu  i  i  A-  t 

,-128:12?) 

t  output  constants  and  the  value  of  pi  . 
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d a  *  a.  :  ow r  t , / z '  000 @083 8 '  / 

dfl  5.  U  6.  6/ 

data  r  ar  add/0, ' 

O  i  /*  i  e  o  r.*  e;  y 

C  -  -  ■  A.  /  A-  _•  • 

c ate  t  *  C-J  t/  6/ ,  1 1  i  n/5/ 


op e r  '  j r.  :  t=2  ,  f  i  1  e— 'mkf  i  1  t .  1  st'  ,  s t at us=/  unkn own  ' ) 

c  Create  output  tile  -for  byte  oriented  I/O. 

cel’  f  z  r  e  a  t  C  '  Enter  output  -file  -for  filter  coef  f  s'  ,  ou  tdes) 


i/.i r  i  te  f 

1 1  ou 

t, 

990 10) 

read  ( 

t  t  :  r, 

j 

99005) 

n  an  g 

wr i t e  1 

ttou 

f 
•  9 

990  28  > 

read  ( 

1 1  i  n 

9 

99005) 

n  term 

wr i  t  e ( 

1 1  o  u 

t  , 

99030 ) 

read  ( 

1 1  i  n 

> 

99006) 

f  whme 

N..I  r  \  m  ( 

♦  tou 

4- 

•  9 

9  y  0  *40  ) 

r  f  3  ( 

1 1  i  n 

9 

99006) 

f  whmb 

I owang  =  /  terms  -  1 >/2 
t  ots.r  -  =  r  a r  g*n  t  e  rm  =• 

990  0^  u--Ti-  f 

9  9  0  g  .£  f  '  IT.  r  t  '  if  2  0  .  1  0  ’ 

99010  *  ; "  t  ;  t  Enter  number  of  angles  '  ,$ ) 

99020  format1''  Enter  number  of  terms  to  compute  each  angle  '  ,* ) 

99830  format <  '  Enter  FWHME  ( ~  6.0)  ' ,%) 

99040  format'"’  Enter  FWHMB  <  ~  1  . 0  >  ',♦> 

■  t  de  =  3 
n  u  m  i  r  t  •  1  0 

n  =  imgdrnt 

:  if  the  number  of  angles  specified  is  too  great  or 

c  n  o  t  e  u  e  r.  1  y  c  i  >•  i  s  :  b  1  e  b  y  four,  exit.  if  the  n  umber  of 

z  angles  is  t oo  gr eat  but  e m e n 1 y  divisible  by  four  this 

o  program,  should  de  modified  by  changing  the  declared  array 
o  sires  to  fit  the  new  requirements.  if  the  number  of  angles 
r  i  s  n  :■  t  e  v  e  n  1  y  d  i  v  i  s  i  b  1  e  by  f  ou  r  this  p  r  ogr  am  w  ill  not 

z  ap  p  1  y  be  z a  u s  e  of  the  extensive  e  x p 1 o i t  a  t i on  of  symme  t r y 

z  r ~  z pert i es  deszr <  bed  above . 

IP  .rang  . gt.  128)  THEN 
wr i te  <  t  tou  t  ,  1  5) 

GOTO  1075 

END  IP 

IP  (mod'.nang,  4)  .  ne.  0)  THEN 

wr i  te  <  t tou t , 20 ) 

GOT  0  10  75 

END  IP 

ndytes  =  nang*4 


co*  up  tad’s.  then  f  nd  the  nurriber  of  rezords  which  will 
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c  be  r  e q  -e; 

r  s  '  ’  -  -  -  -  -f 


!uf  tc  saving  the  number  o-f  pixels  inside  t h ■ 
cp  r  obat  1  i ty  at  each  angle  into  the  lock  up 


IP  *  rr.  c  0  •  <  r  a  r  g*4  ■  ,  256  >  .  e  q  .  O  >  THEN 

r,urr.r  ec  =  numrec  +  ■  nang*4/2561| 

"urrmec  =  numrec  +  ( nang*4/256>  +  1 
END  IF 

c  set  the  an  glee  useful  to  e x p 1 o i ting  the  symme t r y  p r op er  ties. 

angl  =  rang/' 4 
an gisi  —  rang.-  2 

c  -ft  rid  the  standard  deviations  and  their  squares. 

cents.;  =  2 . 0 *sqr  t  (  2 . 0 *a  1  og (  2 . 0  )  > 
s  ;  gm  a  e  =  f  w  h  m  e  /  c  o  n  f  a  c 
s  i  gm  a.  b  =  f  w  h  m  b/  •:  o  n  f  a  c 
s 1 gesq  =  s i gmae*s i gmae 
s  • gb s  q  =  si gmab*s i gmab 

;  -  i  n d  the  r  ; r  eme n  t  f  or  the  f  unc  t  i  on  e  v a  1  u a  t  i  on  dur  i  n g 
c  the  numer i cal  integration  and  the  scale  factor  to  multiply 
;  the  res-'t  by,  wh ’ ch  is  partially  due  to  the  gaussi an  and 
c  partia’ly  due  to  the  integration  method. 

de  1  =  p  i  x/f  I  oa  t  ( n urn  i  n  t  > 

sea' e  =  <de! *del >/<2.0*p i »si gmae*si gmab) 

;  *■  no  the  sem  i  rtiaj  or  and  semi  mi  nor-  axes  of  the  ell  ipse 
;  of  isoprobability  at  the  given  number  of  standard  deviations, 
;  to  be  used  to  del ineate  which  columns  and  rows  to  search 
;  fo"  pixels  inside. 

semaj  =  nostde*s i gmae 
sem  ;  n  =  n os  t  de ^s  i  gmiab 
sernesq  =  semaj  *semaj 
semisq  =  semi  r-*sem!  r, 

■ •!  iTic.X  =  =■  TT*  -s.  .J 

x  m  •  n  =  -  =  •=*  m  r  , 

c  :  1  max  =  :  n  t  <  ■'  xmax/p  i  x  >  +  1,0  ) 
c  ■;  1  m  :  n  =  -  c  o  1  m  a  x 
r  owni  •?.  x  =  c  o  ]  rri  «l  >. 

f  own-  i  r-  =  -  r  oujrn  *.  x 

;  search  for  all  of  the  pixels  inside  the  ell  ipse  at  each 
c  angle  and  tabulate  the  row  and  column  at  which  they  occur, 
c  a! sc,  keep  track  of  the  total  number  for  each  angle. 

tot n urn  =  0 
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wr 
T T>S 


i  tef ttout,  220 1 0  >  an gn um 
t (  '  start  i  n  g  an  g  1  e  ,  i  5 ) 


55y  ,  i  o  :  —  c  o  1  fTi  i  n  ,  col  r rax 
DC  540,  row  =  r owm in,  r owm ax 
: rr, a  t r  x  \  r  ow  ,  col  )  =  6 . 8 
DON'T  I NUE 


Q  ”■  fi  _  : 

r  7  b  TTi  ■ 

t  h  e  t  *  = 


+  «: 


nor  =  -low  an  g,  Iowan  g 

<  p  i  #-f  1  oa  t  <  <  angnum  -  1 )  *n  terms+m  i  nor  )/-f  1  oa  t  (  t  o t  an g)  ) 

< p  i  /2 . 8  > 


a  I  n  t  h  e 
n um i n  a 
DC  108 
xce 
DO 


=  c  o a  ( that  a  > 

=  a  I  n  ( t  h  e  t  a  > 

(angnum)  =  6 
5  col  =  colmin,  col  max 
n  t  =  p i x*f 1  oat ( col  ) 

1018  r ow  =  r  owm in,  r owm ax 

ycen  t  =  -p  i  >: * 1 1  oat  (  r  ow) 

xrot  =  Cxcen t*cos the)  +  ( ycen t*s i n the ) 

yrot  =  C -xcen t *s i n the )  +  < ycen t *c os  the ) 


I F  ( ( x r  o t  . ge .  xm i n )  . an d .  ( x rot 
r  o  t  a  q  =  x  rot*  x  r  o  t 
y  r o  t  a q  =  y  rot*  y  r  o  t 

ymaxsq  =  aem i aq*( 1 . 8  -  ( xr otaq/eemaaq ) ) 
yminsq  =  -ymaxsq 
IF  (Cyrotsq  .  ge  .  yminsq) 


xmax ) >  THEN 


and.  Cyrotsq  .le.  ymaxsq))  THEN 


sisqth  =  s i n the*s i n the 
sincoa  =  a i n the*costKe 
cosqth  =  costhe*costhe 

x  sqrrnj  1  =  -0 . 5*  C  <  a  i  sq  t  h/s  i  gbsq )  +  C  c  osq  t  h/ a  i  ge  a  q  ) ) 
xymu 1  —  C ( si gesq-a i gbsq ) *s i ncos)/< a i gesq*s i gbsq 
ysqmu 1  =  -8 . 5*  C  C  s i sq  t  h/s i ge  sq )  +  C  c  osq th/s i gbsq > ) 
u  =  8.0 

x  1  o  1  i  m=p  i  x  *  <  -f  1  oa  t  ( c  o  1  )  -  0.5) 
y!  ol  im~-p  i  x*C-f  1  oatCrow)  +  8.5) 

DO  1625  xindx  =  1,  numint 

x  =  xlolim  +  <  de  1  *< -f  1  oat  <  x  i  n  dx  >  -  0.5)) 

xsq-fac  =  x  *  x  *  x  sqrriu  1 

DO  1  8  38  y  i  n  dx  =  1  ,  n  um  i  n  t 

>•  =  y  1  o  1  i  m  +  (del  *<f  1  oat  C y  i  n dx )  -8 . 5 >  > 
xytac  =  x*y*xymul 
ysq-f  a;  —  y*y*ysqmul 
u  =  u  +  exp  <  xsq-fac  +  xytac  +  ysq+’ac) 

CONT I NUE 
CONTINUE 

uol p i x  =  u*scal  e 

cms.tr  x(r  ow  ,  c  o  1  )  =  cma  t  r  x  (row,col  )  +  v  o  1  p  i  x 
END  IF 
■ID  IF 
'I  N!JE 
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v o I ume  <  an  gn  um )  =  0.0 
DO  890,  col  =  col  min,  col  max 
DC  880,  row  =  rowmin,  rowmax 

IF  < cma trx(row,col )  . n e .  0.0 )  THEN 
index  =  index  +  1 
r ow tab < i ndex ,  angnum)  =  row 
col tabC i ndex ,  angnum)  =  col 

coe-f-f2(  index,  angnum)  =  cmatrx < row , col  )/-f  1  oat < n terms) 
vol  ume  (  angnum)  =  vol ume ( angnum)  +  coe-ff  2(  i  ndex  , angnum) 
END  IF 


qqg 

CONTINUE 

<5 

890 

CONTINUE 

n um i n  s < an  gn  um )  =  index 

$ 

IF  ( num i ns< angnum)  . gt.  1000)  THEN 

V 

wr i te ( t  tou  t , 25) 

GOTO  1075 

■f. 

thi 

END  IF 

**  • 

c  Use  s 

ymmetry  to  turn  coe-f-ficients  around  a  bit. 

i 

DC  1035  index  =  <(( num i ns< angnum)  +  l)/2)  +  1), 

num i ns( angnum) 

:oe-f-f2<  index,  an  gn  um )  = 

* 

coe-f  -f 2(  <.  i  ndex  -  <  (num i  ns <  angnum)  +  l)/2>) 

,  angnum) 

„•* 

rowtab< i ndex ,  angnum)  =  -r owtab( < i ndex  - 

* 

<< num i ns< angnum)  +  l)/2>),  angnum) 

col tab< i ndex ,  angnum)  =  -col tab( ( i ndex  - 

* 

(< num i ns( angnum)  +  l)/2>>,  angnum) 

1 

1033 

CONTINUE 

totnurn  =  totnum  +  nurri  i  ns <  angnum) 

I  -  <mod<  <.  numi  i  ns<  angnum)  *4)  ,  256)  .eq.  0)  THEN 

r.  umrec  =  numrec  +  <  2*<  num  i  ns(  angnum)  *4/256) ) 
r_qc 

j 

numrec  =  numrec  +  < 2*< < num i ns( angnum) *4/256) 

+  1  )  ) 

r 

END1  I F 

y 

-  :  te < t tou t , 45)  angnum,  num i ns ( angnum) ,  vol ume < angnum) 

1  0  0  0  0 

i  '  d  ce.e  the  number  o-f  contiguous  records  required  to 
:  t  *- e  e-tire  look  up  table. 

-  ■  t e  1 1 o j  t , 30 )  n umr e  c 

c  output  :  -epo-t  indicating  all  o-f  the  input  parameters 
c  selected. 
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num i ns( angnurri)  =  numins(ang2  +  2  -  angnum) 
vol  uP'e  ••  angnum)  =  volume  <ang2  +  2  -  angnum) 

DO  1045  index  =  1,  num i nsC angnum) 

coef f  2(  \  ndex  ,  angnum)  =  coef  f  2(  i  ndex  , 

*  <ang2  +  2  -  angnum)) 

row  tab1!  index,  angnum)  =  col tab< i ndex , 

*  (ang2  +  2  -  angnum)) 

col  tab<  i  ndex  ,  angnurri)  =  row t ab (  i  ndex  , 

*  <ang2  +  2  -  angnum)) 

1645  CONTINUE 

wr  i  te  <  t  tou  t ,  45)  angnurri,  num  i  ns  (  angnum)  ,  vol ume < angnum) 
1040  CONTINUE 

:  -find  the  coefficients,  rows,  and  columns  for  the 
:  last  90  degrees  worth  of  angles  using  the  symmetry  arguments 
:  stated  earlier. 


DO  1050  angnum  =  <ang2  +  2),  nang 

vol ume < angnum)  =  vol ume < nang  +  2  -  angnum) 
num i ns( angnum)  =  num ins( nang  +  2  -  angnum) 

DC  1055  index  =  1,  num i ns< angnum) 

coef f 2 ' . n dex ,  angnum)  *  coef f 2< i ndex , 

*  (rang  +  2  -  angnum)) 

r  ow t ab <  I  nde  >" ,  an gn urn >  =  -rowtab<  i  ndex  , 

*  (nang  x  2  -  angnum)) 

::! tat < i ndex ,  angnum)  =  col tab< i ndex , 

*  (rang  ♦  2  -  angnum,'^ 

1055  CONTINUE 

wr  i  te  ^  t  tou  t  ,45)  angnum,  num*  i  ns(  angnum)  ,  vol  umie  ( angnum) 

[050  CONTINUE 

wr  te  <  1 1 ou  t  ,  50  '•  totn urn 
wr  ;  t  e  <  1 1  ou  t  ,  55  >  n  umin  e  g 


:  save  the  array  containing  the  nurriber  of  pixels  inside 
:  the  ell  ipse  at  each  ang1  e  mto  the  look  up  table.  this  w  i  1  1 
:  be  used  to  mirNmize  the  number  of  evaluations  that  miust  be 
:  performed  in  the  subsequent  program  implementing  the 
:  convo’ ut : on .  note  that  the  number  of  updates  is  not  the 
:  same  for  every  angle. 


t  e ( ou t  de  s ,  n  urn  ins,  n  by  t  e  s  > 


tte  coefficients  into  a.  one-dimensional  arra* 
i tate  transfer  to  disk,  and  calculate  the  address 
t.  then  save  on  disk. 


DO  I960  angnum  =  1 ,  nang 
wr  i  t  e  <  1 1  ou  t ,  90  22  >  an  gn  um 
DO  1065  index  =  1,  num i ns< angnum) 

coeffl< index)  =  coef f2< i ndex ,  angnum) 

c  f  s  e  t  (  i  nde  x  )  =  <  n  *r  ow  t  ab(  index,  angnurri)  )  +  col  tabf  i  ndex  ,a  n  gn  um ) 


021 

f  orma t  ' 

*  . 

i  r  d  =  '  ,  i  6 '  r,  =  ' 
c  *  f  se  t  —  ■"  ,  1  0  ) 

, i 6 '  r  ow  tab  = '  ,  i < 

:  ,  '  col  tab  = 

0  22 
fj  »  = 

-  -  -  -7,  :  *  (  •' 

r  hnt ’ 1 

m  a  k  :n  o  r  :  e  s  . 

-1UE 

s  t  ar  t  .  n  g  an g  1  e 

,  C“  , 

>  •  - 
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r.bl  =  num  \  ns(  angnum)  *4 

cal'!  fwr i te t ou tdes ,  coeffl,  nbl> 

call  fwr  i  te  <  ou  tdes  ,  offset,  nbl  > 

1 6  6G  CONTINUE 

c  -e cord  a  permanent  record  of  the  look  up  table  generation 

c  or  the  pointer. 

write <2, 35)  f whme ,  fwhmb,  pix,  nostde,  nang,  numint,  n 
wr i t  e ( 2 , 30 >  n  umr e  c 
wr i te  < 2 , 48 ) 

DC  1078  angnum  =  1,  nang 

write<2,45)  angnum,  num i ns< angnum) ,  vol ume < angnum) 

1878  C  CN'T I N  U  E 

wr  i  t  e  <  2 , 50 )  t  o  t  n  urn 
wr i  t  e  <  2 , 55  >  n  umne  g 

c  term i mate . 

1075  CONTINUE 

c ’ c  se ( u  n i t=2 ) 


« *. 

A 

re.  1  i  f  c  1  ose  <  1  uou  t 

) 

ta 

stop 

/ 

cr 

format  < f  > 

a 

10 

f  ormat < *  > 

15 

format </ , 

iA 

*  "  number 

of  angles  too  high  for  available  array  space", 

I 

*  "  program 

term i rated . " ,/) 

20 

format </ , 

*  "  rang  miu 

st  be 

evenly  divisible  by  four,  program  aborted.",/) 

£ 

me 

format < " , 

*  "  the  number'  of 

pixels  within  kernel  of  gaussian  error  densit> 

*  "  exceeds 

array 

size,  program  aborted.",/) 

V 

38 

format </ , 

!*• 

*  "  ’ ook-up 

table 

requires", i4,"  contiguous  records.") 

■->  w 

format ( / , 

*  "  longitudinal 

fwhm  =",f7.3,"cm",/, 

£ 

*  "  transverse  fwhm  =  " , f 7 . 3 , " cm" ,/ , 

*  "  pixel  size  =  " , f  6 . 3 , / cm " , / , 

*  "  "  , f 4 . 1 , "  standard  deviations  of  error  density  included", 

*  -■  •  r  !/  ivsl  "  / 

*  "  ’ ook-up  table  generated  for  data  taken  from  ",i3,"  angles.",/, 

*  "  r  i emanr  sum  uses  " , i 3 , "  subd i v i s i ons  per  p  i  xe 1  .  "  ,/  , 

*  "  dimension  of  irriage  array  assumed  to  be  ",i3,/> 
forme  t " ," , 

*  "  angle  number  of  pixels  volume  over  ",/, 

*  "  number  inside  ellipse  these  pixels  "  ,/) 

f  or  ma  t  <  3x ,  i 3 ,  1 3x ,  i 4 ,  1 3x ,  gl 2 . 5) 

forma t < 17x ,  i6) 

format </, i 3, "  negative  coefficients  were  found.",/) 
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*#***¥■*****#*#****#*¥:#*##*#*###*#****#***##*##*#####*#**#*#**#**#* 


PROGRAM  TITLE 


sp  i  nbal  1  .  f 


WRITTEN  BY 
DATE  WRITTEN 
WRITTEN  FOR 


J.  Trent  Woh  1  sch  1  aeger 
3/12/86 

Delay  Doppler  Radar  Imaging 


PROGRAM  INTENT 

""-'1 2  PROGRAM  CREATES  A  N  X  N  REAL  ARRAY 
REPRESENTING  A  SPINNING  BALL.  THE  ARRAY  IS  THEN  STORED  ON  DISK, 
PIXELS  ARE  DIVIDED  INTO  100  <10X10)  SUBPIXELS. 

SUBPIXEL  IS  INCLUDED  IN  THE  REGION,  THEN 
INCREMENTED  BY  1.  THUS  EACH  PIXEL  HAS  A 
0  TO  100. 


IF  THE  CENTER  OF  A 
THE  PIXEL  VALUE  IS 
VALUE  RANGING  FROM 


REG I  ON 


DESCRIPTION 


GEOMETRY 


REGION  1 


SURFACE  OF  CIRCLE 

SPINNING  BALL  CENTER 
RADIUS 


<0 ,0) 

8.00  CM. 


INPUT  FILES 

LU2  =  UNIT  TO  RECEIVE  ANSWERS  FROM 
OUTPUT  FILES 

LIJN  =  IMAGE  FILE  TO  CONTAIN  "  SPINNING-BALL "  TARGET 


£.**  ************  ******  ***#*******#**********#*##*##*#«***•#*#*##***## 


comm,  on  /files/  1  un 
Implicit  Logical  <A-Z> 

REAL  LAM B  DA  <  6 5 5  3  6 > 

INTEGER  N 
i n teger  1 un 
data  7  un  /-!/ 

ca1 1  fcr eat < / Enter  output  file  for  target  model  <sp i n_ball .???)' , 

*  1  un  > 

ppi:-r  *,  ■'  PLEASE  INPUT  THE  IMAGE  ARRAY  SIZE.  <256)  ' 

READ  * ,  N 

FRINT  *,  -  THE  IMAGE  IS  '  ,  N ,  '  BY  '  ,  N,  '  PIXELS . ' 

CALL  b I dbal 1  < 1 ambda  ,  n ) 

END 


ne 


bl dbal 1  (lambda 
t  u  a  1  m a  inline  as 
s  adjustable  ar 


n ) 

a  Subroutine 
ay  dimensions 


to  take  ad v 
for  lambda 


an t  age 


c-c 


c ommon 


t  i  i  e  s/ 
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Implicit  Logical  (A-Z) 

REhL  HAFNP1  j  LAMBDA <  n  ,  n  >  ,  R1 ,  RADSQ ,  R1SQ,  SUM,  X,  XSQ ,  Y,  YSQ ,  PI 
INTEGER  I,  J,  N,  NEYTES ,  i2,  J2,  N2 
integer  1 un 

Parameter  CPI  =  3.14159265) 

DATA  PI  /32 . 8/ 

DATA  R13Q  /I 024.8/ 

HAFNP1  =  8.5  *  FLOAT <N  +  1) 

DO  11  1=1,  N 
DC  15  J  =  1,  N 

LAMBDA < I , J>  =  8.0 
CONTINUE 
CONTINUE 

N2  =  N/2 

DC  1000  1=1,  N2 

PRINT  *,  '  Starting  Column  ' ,  i 
DO  1018  12  =  8 ,  9 

X  =  <1  +  12/10.0  -  HAFNP1 ) 

XSQ  =  X**2 
DO  1020  J  =  1 ,  N 

DC  1038  J2  =  0,  9 

Y  =  -CJ  +  j 2/1 0.0  -  HAFNP1 > 

YSQ  =  Y**2 
RADSQ  =  XSQ  +  YSQ 
IF  < RADSQ  .LT.  R1SQ)  THEN 
PIXEL  CENTER  LIES  WITHIN  REGION  1. 

LAMBDA < I , J)  =  LAMBDA < I , J)  + 

*  <1.1  -  (XSQ  +  YSQ)/R1 SQ) **- 0 . 5*abs(X/Rl ) 

END  I F 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 


1030 
1020 
1010 
1  8  8  0 


set  — HE  T0TAL  INTENSITY. 

SUM  =  0.0 
DC  2020  1=1,  N 
DC  2030  J  =  1 ,  N 

SUM  =  SUM  +  LAMBDA < I , J) 

CONT I NUE 
CONTINUE 

PRINT  *,  "  Total  Intensity  generated:  ",  sum 

SAME  THE  ARRAY  LAMBDA,  WHICH  IS  FILLED  WITH  THE  INTENSITY 
LAMBDA  OF  EACH  PIXEL  FOR  THE  REFLECTION  PROCESS,  ON  THE 
DISK. 


■f  7-  P  (  1  tj  f;  ) 

•  =  =  Ti  *  T:  *  3 
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title 


cwddrre turn  .  f 


Ur  i  1 1  e  n  d  y 
D at*  wr  1 1  e  n 
Wr  i  1 1  e  r,  f or 


J.  Trent  Woh 1 sch 1 aeger 

May  23,  1986 

DDR 


° "  c g r  am  ,  n  ter,  ♦ 

~r,  =  program  convolves  a  set  of  ambiguity  ■functions 
w  th  a  scattering  function  to  produce  a  set  of  deterministic 
De 1 ay-Dopp 1 er  radar  returns,  as  would  be  obtained  by  passing  a 
-adar  signal  through  a  bank  of  BPMF-SLED's. 


:  "  -  *-  files 

lutat  =  look  up  table  containing  the  coefficients  and 
addresses  used  to  implement  the  convolution 
of  a  scatter  function  array  with  the  two-dimensional 
smt  .  g>.  it:/  function.  The  table  has  two  parts. 

The  first  part  is  an  Integer  array  with 
namb  elements,  where  namb  is  the  total 
number  of  ambiguity  functions.  This  array  indicates 
.  how  many  addresses  and  coefficients  are 

present  in  the  second  part  of  the  table  for 
each  ambiguity  function.  The  second  part  of  the  table 
cons' st s  of  the  coefficients  followed  by 
the  addresses,  for  each  ambiguity  function. 

The  table  generator-  is  mkfilt.f. 


♦  -•■■*♦*•***************•*******************•*******#•************«■* 


Imp’ Icit  Logical  <A-Z> 
Include  perm . com" 

I  o-  a-: :  n  file  s/  1  u  n  s  ,  1  u  t  ab 

C'ha-  a:  ter  t-me*24 


: e a ’  s : a 1 1 e -  r t r n s  i  z >  ,  -eturns< imgdmt**2> ,  coeff 15800),  pet 
I'tege-  rcd_es,  am:  bn  urn. ,  colend,  co'st,  i, 

♦  lutab,  namb,  luns  <18),  lutim, 

♦  nsq,  num ; ns< 1 28) ,  of f se t ( 58 00 > , 

♦  p  i  -  end,  pi  x  fir,  pi  x  n  um  ,  p  i  x  s  t  , 

♦  row,  r  owe  n  d  ,  r  ows  t  , 

♦  1 .1,  ;  d  t  hi 


1  r  ‘rCi?'  n b>- 1 e s 

,  ran  add, 

nzero, 

i  or  d 

Integer  1 1  ou  t , 

1 1  i  n 

Integer  length 

, 1 ceoi  ,  s 

t  r- 1  amb  , 

r  t  n  c  n  t 

Data  I 

'  0  0  R  0  0  P  =18 

r  :  -  *•  x  . 

*  *i  t  i  r,  ’  *z 

IA__  +cr  eat  1 '' Enter  output  header  file  name 

2  ALL  fere  a  t  *1  -'Enter  output  Data  file  name 

I-ALL  fopen  (  Enter-  file  containing  target  mode  i 

—  page  43 


'  ,  1  u  n  s  '  8 )  l1 

'  ,  1  u  n  s  1  3 11  > 

■  ,  1  U  "■  r  4  )  . 


IT  7.1 
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i_L  f  rc-enr  ••  "Enter  tile  containing  filters  , lutab) 

:!NT  *,  Enter  number  of  Returns  ' 
ad  (  1 1  :  n  ,  2010)  namb 

r  m  a  t  i  0  '■ 

r t  am t  =  1 

.fit  r  the  *irst  part  of  the  lookup  table,  which 

nc i cates  the  number  of  coefficients  and  addresses  stored 

n  the  second  part  of  the  table,  for  each  ambiguity  function 

r'  g  t  r  -  4*namt 

iLL  fread( lutab,  numins,  length,  lceof,  rtncnt) 

er  < u  n i t= 1 u  t i m , f i  1 e=/ de  t r  t r n . t i m"  , st a tus=' unknown ' ) 

w  i  n  d  1  ■_  t  i  rr, 

iLL  f date^t: me) 

itnurr,  =  1 

t  e  •  1  u  t  :  m  ,  90  0  "  a  - bn  jit.  ,  t  i  me  ( 1  :  24) 


’st  =  1  +  •:  imgdmt  -  imgdim)/2 
■  lend  =  imgdmt  -  (imgdmt  -  imgdim)/2 

;w  St  =  C  o  1  S  t 
i  we n d  =  c o lend 

=  i  mg  dm  t  *  1  mg  dm  t 
t  h  =  - o lend  -  c o  1  s t 
fir  =  >.’  i  mgdmt*<  rowst  -  1  )  )  +  col  st 

lain  loop.  The  cw  algorithm  for  each  Return  is  performed 


r.'TES  =  n  s  n  +  4 
in  add  =  6 

(LL  SYS  I  C  (  PBLK ,  I  ORD  ,  1  u  n  s  <  4  >  ,  scatter,  MBYTES  ,  RANADD  > 

100.  am bn  urn  =  strtamb,  namb 

e  1 1  c  .  t  ,  1  0  ;•  am  bn  urn ,  n  amb ,  t  ime<  12: 19) 

DO  1060  p i xnum  =  1,  imgdmt**2 
re  tu r  n  s ( p i x n um )  =  0.0 
CONTINUE 

length  =  num i ns< ambnum) *4 

CALu  fread(lutab,  coeff,  length,  lceof,  rtncnt) 

CALL  f read' 1 u tab ,  offset,  length,  lceof,  rtncnt) 
wr  i  ter  ft  ou  t  ,  990  0  )  n  urri  i  n  s  (  ambnum ) 

'Tit ' '  non-ze'o  filter  pixels  for  this  Return:  ' , i 5) 
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4.  -  “  =  -  ■£  *  +  J.i.,  ,  [£^  ^  U- 

1  2fr  pin .jfT!  =  p  ■  «  V,  p  i  x e  n c 
IF  scatter1’,  p-i  x  n  um )  .  g  t  .  0.0)  THEN 

r  c  e  '  :  =  n  z  e  r  o  +  1 
DO  1020  i  —  1,  n urn : n s ( ambn um ) 
add' e  s  =  pi  x  n  um  +  of f  se t ( 1  ) 

it  r 
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aid' 


CONTINUE 
END  IF 
CONTINUE 


=  re  turns < addres) 

+  <  s c  a  1 1  e  r  (  p  i  x  n  um  )  * c  oe  i  -f  <  i  ) ) 


pixst  =  p  i  xst  +  imgdmt 

ONTINUE 

ct  =  1  08  .  *+ 1  oat  <  nzer  o)/-f  1  oat  <  i  mgd  i  m**2) 

ir  i  t  e  <  1 1  ou  t ,  ?8  33 )  pet,  i  mgd  i  m  ,  i  mgd  i  m 

ia  t  '.T  x  ,  +'7 . 2 ,  'V.  o-f  pixels  in  the  " , i 3, "  by  -',13, 

"  image  were  nonzero.") 

-f  bate  <  time) 

'P  e  n  <  u  n  i  t= 1  u  t  i  m ,  i  i  1  e= "  de  t  r  t  r  n  .  t  i  mi "  ,  s  t  a  t  u  s= "  u  n  k  n  own " ) 

m  ;  te ( ’ u  t i m , 9007)  ambnum ,  t  i  me  ( 1  :  24) 

iat<"  Begin  processing  Return  #",i2,"  at  ",a> 

I  ose  (  un  i  t=I  u  t  i  m) 

tpu  t  <  re  turns) 


i z  - ose < 1 u tab) 


Return  Construction  algorithm. 

Return  number:  ",i2,"  out  o-f  ",i3,".'  ,/, 

Start  t i me :  ",  a) 
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c  ****************************************************************** 


c  Program  title 

r* 

c  Written  by 

c  Date  written 

c  Written  tor 


rndmr  trn  . f 

J.  Trent  Woh 1 sch 1 aeger 
Nov  5,  1986 

De 1 ay_dopp 1 er  Radar  Imaging  group 


c  Program  intent 

c  This  prograrti  convolves  a  set  of  ambiguity  -functions 

c  with  a  scattering  function  to  produce  a  set  of  random 

c  De 1 ay-Dopp 1 er  radar  returns,  as  would  be  obtained  by  passing  a 

c  radar  signal  through  a  bank  of  BPMF-SLED's. 

c 

c  Input  files 

c  lutab  =  look  up  table  containing  the  coefficients  and 

■:  addresses  used  to  implement  the  convolution 

c  of  an  scatter  function  array  with  the  two-dimensional 

c  ambiguity  function.  The  table  has  two  parts, 

c  The  first  part  i s  an  Integer  array  wi th 

c  namb  elements,  where  namb  is  the  total 

c  number  of  ambiguity  functions.  This  array  indicates 

c  how  many  addresses  and  coefficients  are 

c  present  in  the  second  part  of  the  table  for 

c  each  ambiguity  function.  The  second  part  of  the  table 

c  consists  of  the  coefficients  followed  by 

c  the  addresses,  for  each  ambiguity  function, 

c  The  table  generator  ismkfilt.f. 

c 

c  program  structure 
i  rndmr trn 

c ***************************** ************************************** 

c 

Implicit  Logical  (A-Z) 

Include  'parm.com' 

Common  /files/  1 uns,  1 u  t  ab 

Character  time *24 

Real  scatter(rtrnsiz),  re turnst imgdmt**2> ,  coeff(5000>,  pet 
Complex  ere  turn < i mgdmt**2) ,  b,  gauss 
Integer  irandl,  irand2 

Integer  addres,  ambnum,  colend,  colst,  i, 

*  lutab,  namb,  luns  <10),  lutim, 

*  nsq,  numins(128>,  of f se t ( 5000 ) , 

*  pi x end,  pixfir,  pixnum,  pixst, 

*  row,  rowend,  rowst, 


I n teger 
In  teger 
I n  teger 

Data  It 
Data  t ‘ 
Data  U 


w  i  d  t  h 

n  bytes,  ran  add,  r,  zero, 
1 1  ou  t , 1 1  i  n  » 

1 ength , 1 ceof ,  strtamb, 


i  or  d 


rtnent 


ORD 
t  ou  t 


/x '00000058'/ 
/ 6/ ,  1 1 i n  /5/ 
1 0*-l/ 
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Data  irandl  /322491219/ 

Data  Irand2  /!/ 

1  u  t  i  m  =  1 

CALL  fcreatC' Enter  output  header  file  name  ',luns<8>> 

CALL  fcreatC' Enter  output  Data  -file  name  '  ,luns<3)) 

CALL  fopen  ('Enter  -file  containing  target  model  ',luns(4)) 

CALL  f op enr(' Enter  -file  containing  -filters  ',lutab> 

PRINT  *,  '  Enter  number  of  ambiguity  functions  ' 
read  Cttin,  2010)  namb 
2018  Format ( i 5) 

strtamb  =  1 

c  Read  in  the  first  part  of  the  lookup  table,  which 

c  indicates  the  number  of  coefficients  and  addresses  stored 

c  in  the  second  part  of  the  table,  for  each  ambiguity  function. 

length  =  4*namb 

CALL  fread(lutab,  numins,  length,  1 ceof  ,  rtncnt) 

r 

open  Cun i t=l u  t i m , f i 1 e='rndmt i me . dat ' , status^' unknown ' ) 
rewind  1  u  t i m 
CALL  fdate(time) 
ambnum  =  1 

wr  i  te  <  1  u  t  i  rri ,  9007)  ambnum  ,  t  i me  ( 1  :  24) 
c 1 ose  ( 1 u  t i m) 

colst  =  1  +  (imgdmt  -  imgdim)/2 
colend  =  imgdmt  -  (imgdmt  -  imgdim)/2 
rows  t  =  colst 
r owe nd  =  colend 


nsq  =  imgdmt* imgdmt 

width  =  colend  -  colst 

pixfir  =  ( i mgdmt*( r owst  -  1))  +  colst 

Main  loop.  The  convolution  algorithm  for  each  ambiguity  function 
is  performed  here. 

N BYTES  =  nsq *4 
ran  add  =  0 

CALL  SYSIO  (PBLK,  I ORD ,  luns(4),  scatter,  NBYTES ,  RANADD) 


DO  1005  ambnum  =  strtamb,  namb 

wr  i  te ( t tou t ,  1  5)  ambnum,  namb,  time (12: 19) 
DO  1000  pixnum  =  1,  imgdmt**2 
creturnCp  i  >•: n urn )  =  0.0 
1000  CONTINUE 


length  =  num i ns( ambnum) *4 

CALL  f  re«.d(  1  u  tab  ,  coeff,  length,  keof  ,  rtncnt) 
CALL  freadClutab,  offset,  length,  Iceof,  rtncnt) 
wr i t e ( 1 1 ou  t , 990 0 )  n urn i n s ( ambnum ) 

9900  Formate'  non-zero  filter  pixels  for  this  Return:  ' , i 5 ) 
nzero  =  0 
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pixst  =  pixfir 
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DO  1015  r  ow  =  rowst ,  rowiend 
p i x e  n  d  =  p i x  s t  +  width 
DO  1820  pixnum  =  pixst,  pixend 

Ic  (scatter(pixnum)  .gt.  0.0)  THEN 
nzero  =  nzero  +  1 

CALL  nrmran  (irandl,  irand2,  gauss) 
b  =  scat  ter ( p i xnum) *gauss 
DO  1025  i  =  1,  num i ns < ambnum) 
addres  =  pixnum  +  of-fset(i) 
er  e  turn < addres)  =  ere  turn < addr  es) 

+  <  b*coe-f  i  (  i  ) ) 

CONTINUE 
END  IF 
CONTINUE 

pixst  =  pixst  +  imgdmt 
CONTINUE 

pet  =  1 00 . *f 1  oat < nzero)/f 1  oat < i mgd i m**2) 
wr 1  t  e ( 1 1  ou  t , 90 38 )  pet, i mgd i m , i mgd i m 
For  ma  t  <  1  x  ,  -f  7 . 2  , ' '/.  o-f  pixels  in  the  '  ,  i  3 , "  by  '  ,  i  3  , 

'  image  were  nonzero.') 

CALL  tdate(time) 

op  e  n  <  u  n  i  t=  1  u  t  i  m ,  f  i  1  e= '  r  n  dm  t  i  me  .  da  t '  ,  s  t  a  t  u  s= '  u  n k  n  own ' ) 
wr i te  < 1 u  t i m , 90  07)  ambnum , t ime ( 1 : 24) 

Format <"  Begin  processing  Return  #',i2,'  at  ',a) 
c 1 ose (uni t=l u  t i m) 

DO  3000  pixnum  =  1,  imgdmt**2 

returns  (pixnum)  =  abs  < ere  turn < p i xnum) ) 

CONTINUE 

CALL  output( returns) 

CONTINUE 

CALL  -f  c  1  ose  <1  u  tab) 
s  t  op 


i 


Forma t  < / , 


END 


/ 

5  •  > 


Return  Construction  algorithm. 

Return  number:  ',i2,'  out  o-f  ',i3,'.'  ,/, 


Start  t i me : 
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c  *************************************************  **************** 
c 

c  Program  title  cwasd2d.ftn 


U)r  i  t  ten  by 
Date  written 
Written  -for 


J.  Trent  Woh 1 sch 1 aeger 
Nov .  5 ,  1 986 

De 1 ay-Dopp 1 er  Radar  Imaging 


c 


c 

c 

r 

r 


r 


C 

C 

c 


r 

r 


Prograrri  intent 

This  program  applies  the  conf  i  dence-we  i  gh  t  i  ng 
algorithm  to  the  reconstruction  process  in  De  1  ay-Dopp  1  er 
radar  imaging  to  produce  an  estimate  of  the  desired  image. 

The  "returns"  data,  sorted  by  burst,  and  converted  into 
ordinary  (x,y)  coordinates,  must  be  available. 

Input  -files 

ludat  =  returns  data  in  squeezed  -format 
lutab  =  look  up  table  containing  the  coe-f -f  i  c  i  en  ts  and 
addresses  used  to  implement  the  convolution 
o-f  an  image  array  with  the  two-dimensional 
ambiguity  function.  The  table  has  two  parts. 

The  first  part  is  an  integer *4  array  with 
nburst  elements,  where  nburst  is  the  total 
number  of  bursts.  This  array  indicates 
how  many  addresses  and  coefficients  are 
present  in  the  second  part  of  the  table  for 
each  burst.  The  second  part  of  the  table 
consists  of  the  coefficients  followed  by 
the  addresses,  for  each  burst. 

The  table  generator  is  mkfilt.f. 


prograrri  structure 
cwasd2d 
genr  tr n 


c ******************************************************************* 


common  /rtrnarea/  rtrnpix,  noadd,  noval 
include  "parm.com'' 

character  time *24 

integer  r trnp i x < 1 28) ,  noadd< 10000) 

Real  noval (10000) 

real  r trn<r trns i z ) ,  cwp i mg( i mgdmt**2) ,  coeff<9000) 

integer  addres,  brstnum,  colend,  colst,  i,  lutab,  nburst, 

*  nsq,  numins(l28>,  offset (9000),  pixend,  pixfir,  pixnum,  pixst 

*  rowend,  rowst,  width 

integer  r t  r n  he  ad , r t r  nda t  a , 1 1 ou  t , 1 1 i n 
integer *2  header (128) 

integer  1 u da  t , 1 u  t  ab , 1 u i mg , I  eng t h , i mgf 1 g , 1 c e of ,  s t r  t br  s  t ,  rtncnt 
da  t  a  1 1  ou  t  / 6/ ,  1 1 1 n  /5/ 


c 


I  u  t  i  m 


1 
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? 


ti 
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ii 


J* 
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2000 

2010 


1000 


°0  0  6 


990  0 


CALL  -f  op  enr("  Enter  -file  with  -filters  "  ,lutab) 

CALL  -f  crest  ("  En  ter  output  -file  -for  confidence  weighted  pre  image 
*  1 u i mg) 

CALL  -f  openr  < "  En  ter  returns  header  file  ",rtrnhead> 

CALL  fopenrf "Enter  returns  data  file  " ,rtrndata) 

wr i t  e  < 1 1  ou  t ,  2008) 

read  ( 1 1  i  n  ,  2010)  n  bu  r  s  t 

f ormat < '  Enter  number  of  bursts  'j^) 
format ( i 5 ) 


do  1000  pixnum  =  1,  imgdrrit**2 
cwp i mg( p i xnum)  =  0.0 
con t i nue 


strtbrst  =  1 


Read  in  the  first  part  of  the  lookup  table,  which 
indicates  the  number  of  coefficients  and  addresses  stored 
in  the  second  part  of  the  table,  for  each  burst. 

Do  the  sarnie  for  the  returns  data. 


length  =  4*nburst 

CALL  f r ead< 1 u tab ,  numins,  length,  lceof,  rtncnt) 

CALL  f  reader  trnhead,  rtrnpix,  length,  Iceo-f,  rtncnt) 


open ( un i t=l  u  t  im  ,f  i 1 e=' cwt i me . dat ' , status=' unknown ' ) 
r  ew i n  d  1 u t i m 
CALL  f date (time) 
brstnum  =  1 

wr i t  e  < 1 u  t i m , 90  0 6 )  br s  t num , t i me ( 1 : 24 ) 

formats  Begin  processing  burst  #',i3, 
c 1 ose  < 1 u  t i m ) 


at 


) 


colst  =  1  +  (imgdntt  -  imgdim)/2 
colend  =  imigdrrit  -  < imgdmt  -  imgdim)/2 


rowst  =  colst 
r ow end  =  col  end 


nsq  =  i  mgdrrit*  imgdmt 

w i  d t h  =  col  end  -  colst 

pixfir  =  ( i mgdmt*< rowst  -  1)>  +  colst 


Ma i n  1 oop 
here  . 


The  cw  algorithm  for  each  burst  is  performed 


do  1005  brstnum  =  strtbrst,  nburst 

wr i te < t tou t , 1 5)  brstnum,  nburst,  time (12:19) 

CALL  genr  trn  <  brstnum  ,  rtrndata,  rtrn"1 
wr  i  t  e  (  1 1  ou  t  ,  o c 0  0  >  r, ,jm  i  n  s  (  br  s  t  n  urri ) 
format1''  non -zero  f  M  ter  pixels  for  th-s  burst:  , 

length  =  numi  ns(br  stnum'1  *  4 

CALL  freaddutab,  coeff,  length,  lceof,  rtncnt'1 
CALL  fread(lutab,  offset,  length,  lceof,  rtncnt"1 
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nzero  = 
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0 
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pixst  =  pixfir 
do  1015  row  =  rowst ,  rowend 
pixend  =  pixst  +  width 
do  1020  pixnum  =  pixst,  pixend 

if  <rtrn (pixnum)  .gt.  0.0)  then 
nzero  =  nzero  +  1 
do  1025  i  =  1,  num i ns( brstnum) 
addres  =  pixnum  +  offset(i) 
cwp i mg( addres)  =  cwp i mgC addres) 

*  +  <  r  trn  <  p  i  xnum)  *  coeff(i)) 

continue 

end  if 
continue 

pixst  =  pixst  +  imgdmt 
continue 

pet  =  100 .*f 1 oat(nzero)/f 1 oat( imgdim**2> 
wr i te ( t  tou  t , 9038)  pct,imgdim, imgd i m 
f  ormat  (  1  x  ,  f  7 . 2 ,  'V.  of  pixels  in  the  '  ,  i  3 , "  by  ',i3, 

*  '  image  were  nonzero.') 

update  the  cumulative  sum  for  the  number  of  bursts  so  far. 
In  case  of  system  crash,  we  can  re-start  in  the  middle. 
CALL  frew(luimg) 
length  =  256 

CALL  f wr i te ( 1 u i mg , header ,  1  ength ) 

length  =  4*imgdmt**2 

CALL  f wr i te ( 1 u i mg , cwp i mg , 1 ength ) 

CALL  fdate(time) 

open  < un i t-1 u  t i m , f i 1 e=' cwt i me . dat ' , status=' unknown ' ) 
wr i te<lutim,9007)  brstnum , t i me < 1 : 24) 
formate'  Finished  processing  burst  #',i3,'  at  ' ,a) 
c 1 ose (uni t=l u  t i m) 
continue 

save  the  confidence-weighted  pre image  onto  disk. 

i  mgf  1  g  =  -1 
length  =  256 
CALL  frew(luimg) 

CALL  fwr i te ( 1 u i mg , header , 1 ength ) 

length  =  4* i mgdmt**2 

CALL  fwr i te ( 1 u img, cwp i mg , 1 ength ) 

CALL  f cl ose( 1 u img) 

CALL  f c 1 ose ( 1 u tab) 

CALL  f c 1 ose ( 1 uda t ) 
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formate/ , 

'  CW  reconstruction  algorithm. 
'  burst  number: 


'  i  ' 


•¥> 

* 


.  i  3 . '  out  of 
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* 


Start  time: 


'  ,a) 
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Program  title 


f  i  n  i  mg . f 


Ur i t  ten  by 
Date  written 
Written  for 


J.  Trent  Woh  1  sch  1  aeger- 

Nov .  16,  1 986 

De 1 ay-Dopp 1 er  Radar  Group 


Pr  ogr  ami  1  n  ten  t 

this  program  performs  the  -filtering  -for  the 
reconstruction  -from  a  pre  image  to  a  -final  image. 


type 

type 


0  ==>  confidence-weighted  filtering. 

1  ==>  most-1 i ke 1 y-p os i t i on  filtering. 

2  ==>  convolve  with  a  circularly  symmetric 

gauss i an  to  form  the  desired  image. 


twodim  =  0  ==>  do  not  form  a  two-d i mens i onal  cw 

pre image  by  an  extra  convolution 
with  a  circularly  symme trie 
gauss i an . 

twodim  =  1  ==>  form  a  two-dimensional  cw  pre image 

and  filter. 


The  equation  used  for  reconstruction  is 
d  =  f  *h/g 

h  ere 

d  is  the  four i er  transform  of  the  desired  image, 

f  is  the  fourier  transform  of  the  pre image, 

h  is  the  fourier  transform  of  a  resolution  cell,  and 

g  is  the  fourier  transform  of  the  filter. 

This  program  reads  in  f  (in  the  space  domain), 
takes  the  two-d i mens i onal  fourier  transform,  generates 
the  filter  h/g,  multiplies  f  by  h/g,  and  then 
1  r.  fre.e  fourier  transforms  the  result  to  form  the 
hr-'  mage  ,  which  is  then  stored  on  disc. 

£  type  =  2,  then  g  is  set  to  1,  so  the  effect  is 
tc  cor  vcl'-*e  the  pre  image  (which  should  be  the 
scattering  function  in  this  case)  with  a 
circularly  symme tr i c  gaussian,  to  form  the  desired 
i  IT»5.  » 

p"-  og-  art  structure 
f  i  r,  me 


1  f1  : 29  19S6  finimg.f  Page  2 


cwt  ;1  t ,  exp  bee,  treqx,  freqy,  fwhmb2,  confac, 
fwhmr  ,  g,  h,  hbar  ,  1  p -f  i  ,  ml  f i 1 t ,  twop  i  ,  maxmag,  arg, 
pi,  p i sq ,  r0,  sb2dsq ,  sigbsq,  sigesq,  s i gmab , 
s  i  gmae  ,  s  i  gmar  ,  s  i  grrib2  ,  s i grsq ,  top  i  sq  ,  v  ,  w ,  x.  i  , 
xi  sq,  x  t  1  ,  x  i  2 ,  fwhme  ,  fwhmb ,  r0  ,  p  i  x 

Integer  i,  i  nv  ,  j,  length,  lupimg,  luimg,  Iceof,  rtncnt, 
r,  n m 2 ,  rase os,  twodim,  type,  rtncnt 
Integer  ttin,  ttout 
Integer  *2  header  'll  28) 

EQUIVALENCE  (header<99>,  FWHME),  < header < 1 0 1 ) ,  FWHMB) 

data  pi  /3 . 1 41 59245/ ,  pix  /0.25/ 
data  1 1  i  n  ,'5/  ,  ttout  /  6/ 

-  =  i  m  g  dm  t 

r v  2  =  imgdmt/2 

‘  •'=*  =  2.0*pi 

P  i  sq  =  p i *P  i 

top  i sq  =  2 . 0*p i sq 

confac  =  2 . 6*sqr t < 2 . 0*al og< 2 . 0 ) ) 

twod  i  m  =  0 
f whmb2  =  0.8 

x  i  i  s  e.e 

x 12  =  0.0 
1  p  f  1  =  0.0 


call  -fop  enr  ( "  En  ter  confidence-weighted  pre-image  f  i  1  e "  ,  1  up  i  mg! 
call  f creat< "En ter  output  confidence-weight  file  ", luimg) 


r  0 

= 

P 

i  x/2 . 

0 

wr 

i  te 

\' 

t  tou  t 

,  5) 

r  9 

a  d  ( 

f 

tin, 

10) 

type 

i/ 1  ~ 

f  te 

i. 

ttout 

,15) 

type 

i/.i  r 

i  te 

f 

1 1  o  u  t 

,55) 

r  * 

ad 

ttin, 

35) 

f  whmr 

w  r 

i  t  e 

( 

t  tou  t 

,  60 ) 

fwhmr 

'£•  i 

am  a 

r 

=  fwhrrir/ 

confac 

s i gr  sq  =  s i gmar  *s i gmar 

# 

twodim  *  0 

IF  (type  . e q .  0 )  THEN 
wr i te( ttout ,28  > 
read  (ttin,  10)  twod i m 
wr i  t  e ( 1 1  ou  t , 25 )  twod i m 


I  -  <  twod  irr;  .  eq  .  1  )  THEN 
w  r  i  t  e  (  1 1  o  u  t  ,  6  5 ) 
read  (ttin,  35)  f whmb2 
wr i t  e ( 1 1  ou  t , 70 )  f whmb2 
END  IF 

s i gmb2  =  fwhmb2/conf ac 
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V‘.‘^ 


Dec 


5  17 ;2? 

1936 

finimg.f  Page  3 

Sb2dsq 

=  s  ; 

gmb2*s i gmb2 

wr 1 te< 

t  tou  t 

,75) 

"•■  e  a  d  < 

1 1  i  n  , 

10)  rase os 

i.\  rite  ( 

t  tou  t 

,80)  rase os 

for  mi  at 

f  / 

••  •*■  J 

* 

/ 

please  input  the  type  of  reconstruction  desired. 

y  / 

>'  J 

* 

s 

type  =  0  =>  confidence-weighted. 

* 

/ 

type  =  1  ==>  most-1 i ke 1 y-pos i t i on . 

* 

•* 

type  =  2  ==>  desired  image. 

'  ./» 

* 

/ 

do  not  include  a  decimal  point. 

' 

/ 

' 

* 

/ 

t  yp  e  =  ? 

'  ,/> 

1 1  ( / 


1  0  -format  <  bn  ,  i  3) 
15  -forma 
26  -forma 
* 


i  3? 


f  c  /  , 


is  it  desired  to  -form  a  two-dimensional  confidence-' 


dome,  i  n 


ir  the  pre image  array  arreal  <  ar real 
ar  real  is  c-f  real  t  yp  e  . 


is  in  t h e 


ace 


/ 


* 

weighted  array  before  f i 1 t e r i n g? 

■/ 

twodim  =  0  ==  >  do  not  forrri 

the  2d  cw 

array . 

'  ,/, 

* 

/ 

twodim  =  1  = >  form  the  2d 

cw  array. 

* 

/ 

do  not  include  a  decimal  point. 

* 

/ 

'  >/» 

* 

y 

twodim  =  ? 

'  ,/) 

25  f 0- mat 

J 

twod  i  mi  =  '  ,  i  3) 

• 

35  forrriat 

<  f  20 

.10) 

55  f o'- rti at 

/  .«• 

•  ? 

* 

please  input  the  full  width  half 

max  form 

of  the 

*  / 

-4- 

/ 

standard  deviation  of  the  circularly  symmetric 

'  |/, 

* 

y 

gaussian  resolution  cell,  f whmr . 

include 

the 

' 

* 

/ 

de  c i mal  point. 

'  ,  / , 

Hr 

/ 

'  ./, 

* 

fwhmr  =  ?  cmi. 

'  »/> 

6  *3  r  3  r  iTi  *.  *i 

fwhmr  —  fl0.5) 

if  -format 

5 

*■ 

/ 

please  input  the  full  width  half 

max  form 

of  the 

/  / 

*■ 

standard  deviation  of  the  circularly  symmetric 

'  ,  / , 

gaussian  which  is  to  be  used  to 

create  the 

two- 

/ 

♦ 

dime n s i on a  1  confidence -we i gh  t  e  d 

pre image  array. 

'  >/ > 

♦ 

include  the  decimal  point. 

/  / 

♦ 

' 

f  whmb2  =  ?  cmi . 

'  ,/> 

fjrrriit 

<  /  ' 

•  > 

f whmb2  «=  '  ,  f  1  0 . 5) 

75  ^.rrr^t 

<  /  » 

* 

/ 

is  it  desired  to  use  a  low  pass 

filter  of 

the 

' 

* 

✓ 

raised  cosine  type? 

/  / 

5  '  5 

* 

/ 

rase  os  =  0  ==>  do  not  use 

this  lpf. 

y  / 

j  y 

*■ 

/ 

rase os  =  1  ==>  use  this  lpf  type. 

> /  * 

♦ 

y 

d:  not  : no’ ude  a  decimal  pci n  t  . 

>  « 

♦ 

r  a  s  c  os  =  ? 

?  •  5 

80  ^orr:s* 

rs.scos  =  ,  i  3)  Cr 
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i-T* 


/V  Vi 
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c«  ■  '  ew(  ’  up  i mg) 

i  .  _  _  4-  w  _.  -“}  cr  y 

§r^tii  —  ^  j  o 

I  P  ( type  .  ne  .  2>  call  f  readU  up  i  mg,  header  ,  1  ength  ,  1  ceof  ,  r  tncn  t ) 
'  eng*!":  =  4*  i  mgdmt**2 

cal"  -f  r  e  ad  (  1  u p  i  mg ,  ar  r  e  a  1  , 1  en  g  t  h  ,  1  c e of ,  r 1  nc n  t  > 

Ic(lcsof.eq.l)  THEN 

wr i te ( t tou t , *>  "  EOF  encountered  while  reading  input  data 
END  IF 

PRINT  *,  'Enter  FWHME' 

READ  * ,  f  whme 
PRINT  *,  "Enter-  FWHMB" 

READ  *  ,  f  whmb 
sigrmae  =  f whme/conf ac 
s  i  gmab  —  fwhmb/conf  ac 
=■  i  ge  sq  =  =■  i  gmae*s  i  gmae 
=  i  gt  a  ;  =  £■  i  gmab*s  i  gmab 

1“  (  r  ascos  .  eq  .  1  )  THEN 

•  ' 2  =  (si gmab/s i gmar ) / < 2*p i x  > 
x  i  1  =  6 . 8  *  x  i  2 
END  IF 

DC  990  i  =  1 ,  imgdmt 
DO  995  j  =  1 ,  imgdmt 
ar i mag( i , j )  =  0 . 0 
CONTINUE 
CONTINUE 

take  the  two-d  i  mens  i  onal  fourier  trans-form  of  the  pre  image. 


a":  *+ 1  <  arreal  ,  arimag,  n,  n,  1,  inv) 
a1'  f -^t  (arreal,  arimag,  1,  n,  n ,  inv) 

:  a  .  x  a  g  =  0 

'C  C^C  =  1  ,  imgdmt 
I  C  271  j  =  I  ,  i  mg  dm' t 

maxima  g  =  max  (rr.axmag ,  sqr  t  (  arreal  (  i  ,  j  )  **2  +  arimag(i,j>* 


~  -  *_■  —  1  ,  '  m  gdm  t 

DC  331  j  =  1,  irrigdmt 


j ’ r  t  a  ( i  , J )  *  arreal  ( i  , j ) /maxmag 
5- ■ mag ' i , j )  =  ar i mag  ( i , j ) /maxmag 
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1010  =  1  ,  mgdm  t 

F  1  ■  .  1  e  .  r.v2  ■  "HEN 

f  r  eq  ■.  =  f  ’  oa  t 1  -  1  )/(  f  1  oa  t  ( n  >  *p  i  x  > 


jViVt 


i  tint  mg . f  Page  5 

=  + 1  oa t  <  i  -  1  > /" ( 4 1  oa t  <  n  )  *p  i  x  )  -  1.0/p  i  x 

j  =  1,  imgdmt 
.  1  e  .  nv2>  THEN 

-f  r  e  q  y  =  4  1  oa  t  ( j  -  1 )  /  <  -f  I  oa  t  ( n )  *p  i  x  ) 
e  1  =  e 

■freqy  =  +loat(j  -  1  >/< -H  oat  <  n  >  *p  i  x )  -  1.0/pix 
END  IF 

x  i  sq  =  <  4 reqx*fr eqx  )  +  ( -f  r  eqy*-f  r  eqy  > 

>■  i  =  sqr  t  ( x  i  sq ) 

IF  (type  .eq.  0)  THEN 

v  =  top i sq*( s i gesq  -  sigbsq)*xisq 
cw-f  i  1  t  =  exp  (  -top  i  sq*x  i  sq*(  s  i  grsq 

-  2 . 6 *s i gbsq ) > /expbe s< v ) 
ar r  e a  1  (  i  ,  j  )  =  arrea  1  <  i  ,  j  )  *c w-f  i  1  t 
a'  mag ( i  , j  >  =  ar i mag < i  , J  ) *cwf i 1 t 
END  IF 

IF  (type  .eq.  1)  THEN 

h  =  exp ( -top i sq*s i gr sq*x i sq > 
v  =  pisq*(sigesq  -  sigbsq)*xisq 
w  =  twop i *x i *r0 
IF  ( x i  .ne.  0.0)  THEN 

g  =  exp < -top i sq»s i gbsq*x i sq) 

*  expbes(v)*besj 1 (w)*r0/x i 

e  ‘  s  e 

g  =  exp < -top i sq»s i gbsq*x i sq ) 

*expbes<  v  )  *p  i  *r0*r0 


END  IF 
m }  4  i  1  t  = 

P 

i *r@*r0*h/g 

arr  ea’ (  i 

}  j 

>  =  arrea 1 ( 

ar i mag< : 

i 

5  v 

>  =  arimag( 

END  IF 

IF  =type  .eq.  2)  THEN 

ar g  =  -top i sq*s i or sq*x i sq 
Ir  (arg  . ge .  -70 T@)  THEN 
h  =  exp (arg) 

e  I  if 

h  =  6.0 

END  IF 

ar  r  e  a  1  (  i  ,  j  >  =  ar  r  e a  1  (  i  ,  j  )  *h 
ar  imiag(  i  ,  j  )  =  ar  i  m a g <  i  ,  j  >*h 
END  IF 

IF  twpdin-  .eq.  1)  THEN 

Lty  =  e:.p  ' -t  Dp  sq*sb2dsq*  -  i  sq  ^ 
arrea  I  ;  , ,  >  =  ar  r  e  a.  I  (  i  ,  ,1  s  *  h  fc  ar 

=*. "  .  rr i  s».  5 1  i  |  v  1  =  s.  r  »  mag' i , J ) *htar 
END  IF 

IF  'rase os  .eq.  1)  THEN 
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END  IF 

IP  <  <  X ■  . ge .  X i 1 >  .and.  <xi  .It.  xi2>)  THEN 

1  p  -f  1  =  0 . 5*  <  1  +  c os  ( p  i  *  <  x  i  -  x  i  1  )  /  ( x  i  2  -  x  i 
END  I  F 

Ic  (xi  . ge .  x i 2)  THEN 
Ipfl  =  0.0 
END  IF 

ar  r  a  a  1  >'  i  ,  j  )  =  ar  r  e  a  1  <  i  ,  j  )  * 1  p  f 1 
ar  i  magC  i  ,  j  >  =  ar  imag(  i  ,  j  )*1  pf  1 
END  IF 
CONTINUE 
CONTINUE 

real  and  arimag  are  now  the  real  and  imaginary  parts  of 
e  Fourier  transform  of  the  desired  image. 

r?  the  inverse  two-dimensional  f^urier  transform  to 
tar  the  desired  image. 


I  ffttarreal,  arimag,  n,  n,  1,  inv) 

1  f -ft  tar  real  ,  arimag,  1,  n,  n,  inv) 

"eal  is  now  the  desired  image  in  the  space  domain, 
ve  the  desired  image  on  disc. 

frew(luimg) 

;th  =  2-- 

;  fwr  :  te <  1  u i mg  ,  header,  length) 
gth  =  4*imgdmt**2 
1  fwr i te ( 1 u i mg , arreal , 1 ength ) 

1  f  c  1  ose  (.  1  up  i  mg) 

1  f  c  1  ose  (.  1  u  img) 

p  '  final  image  written  to  disc  ' 
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r  #**■***#***#*#*****#***#**#######■##############*######*#*####**##### 


PRO 3 RAM  TITLE 
WRITTEN  by 


u  u  f  i  n  1  .  f 


'EN 


J.  Trent  Woh 1 sch 1 aeger 
Nov  18,  1986 

De 1 ay-Dopp 1 er  Radar  Imaging  Group 


’  R L i RAK't  I  NT ENT 

-Ms  is  a  graphics  utility  program  which  displays 
-final  image  created  by  the  cw  algorithm. 


r  **#•*•*•»**•**  *##*##***####*#**#***###*****#*#****#***#***#***#* 


implicit  Logical  <A-Z> 

Include  '  p arm . c om ' 

Integer  xof  f  ,  yof  f  ,  xs  i  z  ,  ys  i  z 
Rea.1  anneal  (imgdmt,  imgdmt) 

Integer  length,  luimg,  Iceo-f,  rtncnt 
Integer- *2  header- (128) 

Integer-  i  ,  j 
Rea’  pic  (128,128) 


C  PARAMETERS 

Parameter  (xsiz  =  0,  ysiz  =  0) 


C  INITIALIZE  GRAPHICS  PROCESSOR 
CALL  M3IASNGP  (0,0) 

C  ENABLE  ALL  PLANES 

CALL  MG'IPLN  (-1) 

CALL  fopenr  ("Enter  conf  i  dence-we  i  ghted  image  file",  luimg) 

CALL  fr-ew  (luimg) 
length  =  256 

CALL  tread  (luimg,  header,  length,  Ice of,  rtncnt) 
length  =  4*imgdmt**2 

CALL  tread  (luimg,  arreal ,  length,  Iceof,  rtncnt) 

IF  (Iceof  . eq.  1)  THEN 

=,RINT  t ,  y  End-of-file  encountered  while  reading  input  data. 


20 
1  0 


pic  (i,j)  =  arreal  (i+64,  j+64) 
CONTINUE 
CONTINUE 


-  :  -  i  t-Q  i 

.  r  '  -  J  *  —  W  f  I 


o-f  i  |  y  o  *f  ■f  j  x  s  i 


- 


